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This  dissertation  presents  a fundamental  study  of 
stress  fields  around  an  elliptical  hole,  with  and  without 
cracks  emanating  from  it, in  unidirectional  composite 
materials.  The  first  part  deals  with  the  stress  concentra~ 
tion  factors.  The  second  part  discusses  the  stress  intensity 
factors  and  energy  release  rate.  Analytical  study  in  the 
first  part  is  based  on  the  closed  form  of  Savin's  exact 
solution  and  the  finite  element  solution  (displacement 
model) . Solutions  of  different  anisotropic  materials  with 
^iffsrent  elliptical  holes  in  shape  and  size  are  given. 

The  fr'acture  mechanics  part  starts  with  a discussion 
of  the  application  of  different  energy  principles  to  the 
finite  element  technique.  This  discussion  is  followed  by  the 


xii 


introduction  of  a new  mathematical  formulation  of  a finite 
element  model.  This  model  is  based  on  the  modified  comple- 
mentary energy  principle  in  conjunction  with  the  anisotropic 
elasticity  and  complex  variable  techniques. 

Calculations  of  stress  intensity  factors  for  single 
mode  and  mixed  mode  fractures  are  included.  Analytical  and 
experimental  evaluations  of  the  strain  energy  release  rate 
are  determined.  Analytical  evaluation  of  the  energy  release 
rate  is  based  on  the  finite  element  solution.  The  experi- 
mental part  of  this  evaluation  is  based  on  the  compliance 
calibration  technique. 

The  finite  element  computer  code  developed  here  seems 
to  offer  accurate  results  in  exceptionally  efficient  time. 
This  code  is  capable  of  handling  the  two-dimensional  fracture 
mechanics  problem  for  isotropic  and  anisotropic  cases  with 
any  niimber  of  cracks. 
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CHAPTER  1 
INTRODUCTION 

1.1  Motivation  of  the  Research 

The  field  of  composite  materials  has  experienced 
tremendous  growth  since  the  introduction  of  the  so-called 
"advanced  composites"  in  the  early  '60s.  This  growth  has 
largely  been  the  result  of  a desire  to  apply  these  high- 
strength  and  high-modulus  but  lightweight  materials  in  air 
transportation  systems  and  national  defense  as  well  as  the 
automotive  industry.  This  weight  reduction  factor  should 
result  in  better  energy  conservation. 

While  the  investigation  and  prediction  of  the 
elastic  fracture  behaviors  of  conventional  engineering 
materials  are  fairly  well  documented  and  understood,  the 
prediction  of  the  fracture  characteristics  of  anisotropic 
materials  have  presented  engineers  with  some  challenging 
problems . 

Prediction  of  failure  due  to  the  presence  of  small 
cracks  or  flaws,  especially  in  composite  materials  containing 
holes,  has  proven  to  be  particularly  difficult.  Therefore, 
in  order  to  overcome  this  difficulty  we  have  to  provide  some 
reliable  tools  to  serve  engineering  needs  and  to  assess  the 
safety  of  structural  components.  In  making  predictions  of 
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this  failure,  we  can  facilitate  the  optimization  of  com- 
posite materials  in  terms  of  fracture  mechanics  parameters 
by  estimating  the  high  stresses  near  the  crack  tip. 

Elastic  fracture  mechanics  provides  such  tools  to  deal 
with  these  phenomena.  A characterization  of  the  crack  with 
its  local  stress  field  by  means  of  the  stress  intensity 
factor  Kj,  and  for  different  loading  modes  (to  be 

discussed  in  the  next  section)  can  fulfill  these  needs. 

Cracks  emanating  from  holes  represent  one  of  the 
most  common  sources  for  catastropic  failure.  According  to 
Karlsson  and  Backland  [1],  approximately  one-third  of  all 
cracks  were  usually  generated  by  high  stress  concentrations 
near  the  cutouts.  In  unidirectional  composite  materials,  the 
location  of  these  high  stress  concentrations  not  only  depends 
on  the  type  of  load  and  the  geometry  of  the  hole,  but  also  on 
the  fiber  orientation  with  respect  to  the  geometric  axes. 

When  attempting  to  analyze  such  a problem  (especially 
with  cracks)  using  an  analytical  method  such  as  the  conformal 
mapping  technique  to  obtain  an  exact  solution,  we  found  that 
it  was  difficult  to  achieve.  In  more  and  more  engineering 
situations  today,  it  is  necessary  to  obtain  approximate 

numerical  solutions  to  problems  rather  than  exact  closed-form 
solutions . 

The  finite  element  method  is  a powerful  numerical 
analysis  technique  for  obtaining  approximate  solutions  to  a 
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wide  variety  of  engineering  problems.  This  technique  was 
first  established  in  the  early  '60s.  Since  that  time  it  has 
gained  more  popularity  and  recognition  and  has  become  an 
acceptable  engineering  analysis  tool. 

The  desire  to  contribute  to  the  further  development 
of  the  finite  element  method  was  one  motivation  of  this 
research.  The  necessity  of  analyzing  the  fracture  mech- 
anics problem  in  unidirectional  reinforced  composite  material 
was  the  other  motivation. 

1.2  Fracture  Mechanics:  A Background  Review 

The  field  of  fracture  mechanics  began  in  1921  with 
the  introduction  of  Griffith's  theory  [2].  This  theory 
states  that  a crack  will  begin  to  propagate  if  the  elastic 
energy  released  by  its  growth  is  greater  than  the  energy 
required  to  create  the  fracture  surfaces. 

In  1958  Irwin  extended  Griffith's  theory  [2]  into 
linear  elastic  fracture  mechanics.  He  pointed  out  that  there 
exist  three  kinematically  admissible  crack  extension  modes 
(Figure  1.1).  The  first  mode  is  characterized  by  the  motions 
of  the  crack  surfaces  that  tend  to  separate  symmetrically 
with  respect  to  the  plane  occupied  by  the  crack  prior  to 
deformation.  It  is  called  the  opening  mode  (Mode  I).  The 
second  mode  concerns  local  deformation  in  which  the  crack 
surfaces  glide  over  one  another  in  opposite  directions  but 
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Mode  I 

(Opening  Mode) 


Mode  II 
(Sliding  Mode) 


Mode  III 
(Tearing  Mode) 


Figure  1.1:  Linear  elastic  fracture  mechanics  modes 

in  the  same  plane.  It  is  called  the  sliding  mode  (Mode  II). 
The  final  mode  refers  to  the  antiplane  shear  in  which  the 
movement  of  the  crack  surfaces  can  be  related  to  the  warping 
action  of  noncircular  cylinders  under  torsion  where  the 
initially  material  points  in  one  plane  occupy  different 
planes  after  deformation.  It  is  called  the  tearing  mode 
(Mode  III) . Each  of  these  fracture  modes  has  been  mathe- 
matically modeled  and  well  defined  in  terms  of  an  important 
fracture  parameter  called  the  stress  intensity  factor  (K) . 
These  stress  intensity  factors  are  the  measure  of  the 
hypothetical  elastic  stress  singularity  at  the  crack  tip, 
according  to  the  classical  theory  of  elasticity. 
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It  is  of  considerable  importance  in  many  engineering 
problems  to  know  the  numerical  value  of  the  stress  intensity 
factors,  since  there  are  critical  values  of  K which  deter- 
mine whether  or  not  the  crack  will  propagate.  In  composite 
materials,  sometimes  the  crack  propagates  in  a different 
direction  from  the  original  crack  direction.  This  leads  us 

to  the  following  brief  discussion  of  some  important  fracture 
criteria. 

The  first  approach  is  that  of  Sih  and  Chen  [3]  in 
which  they  formulate  a strain  energy  density  function,  S,  a 
measurement  of  the  elastic  energy  field  in  the  vicinity  of  a 
crack  tip.  This  function,  S,  is  a quadratic  form  of  the 
Mode  I and  Mode  II  stress  intensity  factors  and  varies  with 
the  angle,  6,  measured  from  the  plane  of  the  crack.  It  is 
then  postulated  that 

1.  Crack  initiation  will  occur  at  the  crack 
tip  in  a radial  direction  along  which  the 
strain  energy  density,  S,  is  at  a minimum. 

2.  The  crack  will  begin  to  propagate  when  the 
density,  S,  reaches  some  critical  value, 

Scr»  which  is  considered  to  be  a material 
constant  determined  from  experiment. 

The  second  approach  of  fracture  criteria  is  postu- 
lated as  follows  [4] : 
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1.  The  crack  will  propagate  in  the  direction 
which  causes  the  maximum  energy  release 
rate  to  occur. 

2.  The  fracture  will  initiate  when  the  release 
rate  in  that  direction  reaches  a certain 
specified  level. 

The  other  approach  to  fracture  criteria  is  to  employ 
the  method  of  compliance  calibration  to  determine  the  strain 
energy  release  reate  experimentally.  This  method  will  be 
discussed  further  in  Chapter  6. 

The  critical  value  of  the  stress  intensity  factor  is 
considered  to  be  the  fracture  toughness.  In  spite  of  being 
considered  as  a material  property  in  most  of  the  conventional 
engineering  materials,  this  fracture  toughness  is  still  a 
questionable  property  in  composite  materials. 

Analysis  of  stress  intensity  factors  at  the  crack 
tip  has  been  the  subject  of  many  numerical  applications. 

Some  of  these  numerical  schemes  are  methods  including  boundary 
collocation  [5-7],  conformal  mapping  [8],  the  analytical  use 
of  Green's  function  [9],  and  the  finite  element  method.  A 
literature  review  in  the  applications  of  the  finite  element 
method  in  fracture  mechanics  is  introduced  by  Gallagher  [10] . 

1.3  Objectives  and  Assumptions 

The  major  goal  of  this  research  was  to  obtain  a solu- 
tion for  the  stresses  near  the  crack  tip  for  the  problem  of 
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cracks  emanating  from  a hole  in  unidirectional  fiber- 
reinforced  composite  materials.  The  secondary  goal  was  that 
of  developing  a special  finite  element  computer  program  to 
treat  this  crack  tip  singularity  in  an  accurate  and  effi- 
cient way.  The  experimental  verification  for  a simple  case 
was  the  final  goal. 

The  following  considerations  are  the  main  basic 
assumptions  throughout  this  research: 

1.  The  composite  materials  are  viewed  from  a 
macroscopic  scale,  i.e.,  to  be  characterized 

^ as  homogeneous  and  rectilinearly  anisotropic. 

2.  The  stress  singularities  associated  with  the 
cracks  are  of  the  order  r , where  r is 
the  distance  measured  from  the  crack  tip. 

3.  There  are  no  thermal  stresses  or  d3niamic 
loading  effects  on  the  crack  initiation. 

4.  There  are  small  elastic  strains  during  the 
deformation . 

5.  There  is  two-dimensional  plane  stress 
elasticity. 

1.4  A Brief  Review  of  Related  Literatures 
The  increased  demand  for  determining  accurate  pre- 
dictions of  stress  fields  in  complex  structures,  such  as  the 
problem  of  cracks  emanating  from  holes,  has  been  the  primary 
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concern  for  many  investigators.  This  section  will  discuss 
some  of  the  work  related  to  that  type  of  problem. 

In  1956  Bowie  [8]  started  one  of  the  first  analyses 
of  this  problem  by  solving  the  problem  of  equal  cracks 
amanating  from  a circular  hole  in  an  isotropic  infinite 
plate.  The  plate  in  his  problem  was  subjected  to  uniaxial 
or  biaxial  stresses  at  infinity  while  the  hole  and  crack 
boundaries  were  assumed  to  be  free  of  tractions.  His  analy- 
sis was  based  on  the  conformal  mapping  of  the  inner  star- 
cracked  boundary  onto  a simple  circle.  Unfortunately,  such 
a mapping  for  a realistic  case  of  finite  plate  geometries 
would  be  extremely  difficult  to  find.  As  cited  in  Oladimeji 
[11],  the  shortcoming  of  this  method  led  Bowie  in  1965  to 
apply  a mapping-collocation  technique  to  solve  this  problem 
for  a finite  sheet  under  uniaxial  loading  condition. 

In  1971  Bowie  and  Freese  [12]  investigated  the  problem 
of  a central  crack  in  a rectangular  sheet  of  orthotropic 
materials  loaded  under  uniaxial  tension.  They  extended  the 
method  of  modified  mapping-collocation,  which  was  originally 

formulated  for  plane  isotropic  analysis,  to  solve  the  problem 
of  plane  anisotropy. 

In  1971  Gandhi  [7]  treated  the  more  general  case  of 
arbitrary  orientations  of  the  crack  and  the  material  principal 
symmetry  axes  relative  to  loading  direction.  He  used  the 
modified  mapping-collocation  technique  to  calculate  the  stress 
intensity  factors  for  several  orientations  of  cracks. 


9 


In  1974  Hsu  [13]  obtained  the  stress  intensity  fac- 
tors at  the  root  of  a radial  crack  emanating  from  a circular 
hole  in  an  infinite  sheet,  made  of  isotropic  material,  under 
inclined  uniform  tension.  He  used  the  Muskhelishvili  formu- 
lation and  the  conformal  mapping  method  in  conjunction  with 
the  truncation  technique  to  yield  results  for  the  combined 
mode  fracture,  and 

In  1980  Tirosh  [14]  used  the  energy  release  rate  to 
predict  the  fracture  strength  of  the  central  crack  in 
t*^i*^i^®ctional  fiber-reinforced  composite  materials.  He 
worked  the  problem  of  a crack  in  the  fiber  direction  with 
the  fiber  direction  being  inclined  to  the  loading  direction. 
He  then  used  mixed  mode  fracture  analysis  to  solve  for  Kj 
and  Kjj. 

In  1980  Wang  and  Yau  [15]  introduced  an  analytical 
study  of  cracks  emanating  from  a circular  hole  in  an  off- 
axis  unidirectional  fiber-reinforced  composite.  They  used 
the  well-known  path- independent  J-integral  for  an  arbitrary 
path  which  begins  on  one  crack  surface,  encloses  the  crack 
tip,  and  terminates  on  the  opposite  surface  of  the  crack. 

Dissertation  Overview  and  Contributions 

The  following  areas  have  been  investigated  in  detail 
throughout  this  work: 

1.  Review  of  Savin's  classical  solution. 
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2.  Review  of  some  variational  principles  and 

their  applications  to  the  finite  element 
method, 

3.  Development  of  a new  mathematical  formula- 
tion in  application  of  the  finite  element 
method  to  fracture  mechanics. 

4.  An  experimental  verification  of  a simple 
case . 

The  first  area  covered  in  this  work  is  a review  of 
the  classical  elliptical  problem  (an  elliptical  hole  in  the 
center  of  an  anisotropic  infinite  plate  loaded  under  the 
state  of  simple  tensile  stress)  which  was  solved  in  closed- 
form  solution  by  Savin  [16]  in  1939.  This  solution  is 
presented  in  Chapter  2.  The  objectives  of  this  review  are 
to  achieve  a full  understanding  of  this  analytical  technique, 
which  can  lead  to  a similar  formation  in  adapting  a new 
numerical  scheme.  In  addition,  Savin’s  closed-form  solution 
has  been  used  for  checking  the  numerical  method  for  the 
problem  of  an  elliptical  hole  in  anisotropic  media.  The 
other  reason  for  presenting  this  classical  solution  is  to 
examine  and  compare  the  stress  intensity  factor  in  different 
composite  materials  under  the  same  loading  condition. 

The  second  area  covered  is  a review  of  some  existing 
and  available  methods  in  variational  principles  (Chapter  3) 
and  a discussion  of  their  application  in  finite  element 
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analysis.  A review  of  the  literature  concerning  such  appli- 
cations indicates  that  most  of  the  work  in  this  area  falls 
into  two  main  categories.  One  is  for  the  characterization 
of  fracture  mechanics  behavior  in  isotropic  materials,  where 
the  presence  of  a single  and/or  a mixed  mode  is  included. 

The  other  category  involves  the  study  of  fracture  mechanics 
in  composite  materials.  In  this  category,  the  literature 
has  emphasized  a single  mode  loading  and  the  case  of  ortho- 
tropic materials,  when  the  material  principal  axes  coincide 
with  the  geometrical  axes,  with  the  crack  in  one  of  the 
material  principal  directions.  By  comparison,  mixed  mode 
fracture  and  generally  anisotropic  behavior  have  received 
little  attention. 

In  this  work  a new  modified  variational  principle 
has  been  extended  and  presented  in  Chapter  4.  This  technique 
is  then  used  to  solve  the  problem  of  cracks  emanating  from 
a hole  in  a unidirectional  fiber-reinforced  composite.  This 
mathematical  formulation  and  finite  element  computer  code, 
which  is  listed  in  Appendix  E,  are  the  main  contributions 
of  this  work. 

In  Chapter  5 we  present  some  results  from  Savin's 
solution  (exact  solution)  and  the  finite  element  solution 
(approximate  solution)  for  some  practical  problems. 

In  Chapter  6,  an  experimental  program  carried  out 
to  seek  some  comparison  between  the  numerical  prediction  of 
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the  energy  release  rate  and  its  experimental  evaluation  for  a 
simple  case  is  discussed. 

Finally,  a discussion  and  conclusion  is  given  in 
Chapter  7,  m addition  to  some  suggestions  for  further  work 
and  research  related  to ‘this  area. 


CHAPTER  2 

REVIEW  OF  CLASSICAL  SOLUTION 
BY  USING  THE  METHOD  OF  COMPLEX  VARIABLES 


The  main  purpose  of  this  work  is  to  study  the  problem 
of  cracks  emanating  from  an  elliptical  hole  in  a unidirec- 
tional fiber-reinforced  composite  material.  We  start  from 
the  general  elasticity  formulation  of  a related  but  simple 
problem.  The  problem  is  to  determine  the  stresses  in  a 
homogeneous  anisotropic  plate  weakened  by  an  elliptical  hole 
at  the  center  of  the  plate  subjected  to  the  application  of  a 
uniform  tension  stress  of  intensity  P at  infinity. 

The  solution  of  this  problem  will  be  used  as  the 
datum  line  for  checking  the  numerical  program,  when  we  apply 
it  to  such  a simple  case.  The  numerical  program  used  in 
this  research  is  to  apply  general  elasticity  theory  in  con- 
junction with  the  principle  of  minimum  potential  energy  and 
the  modified  principle  of  minimum  complementary  energy  to 
formulate  a finite  element  technique. 

The  unidirectional  composite  material  is  considered 
as  an  orthotropic  material.  In  general  we  have  six  inde- 
pendent constants  when  we  apply  Hook's  law  in  a three- 
dimensional  element.  There  are  four  independent  constants 
for  a two-dimensional  problem  which  will  be  our  case 
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of  study.  For  the  case  when  the  fiber  direction  is  not 
parallel  to  geometric  axes,  there  are  six  elastic  constants, 
which  can  be  expressed  in  terms  of  the  original  four 
constants.  Therefore,  for  the  two-dimensional  case  there 
are  only  four  independent  material  constants. 

In  this  chapter  the  formulation  of  the  elasticity 
problem  using  Savin's  procedure  [17]  for  the  complex 
variables  representation  of  the  classical  theory  of  elas- 
ticity will  be  discussed.  Then  the  solution  for  the  stress 
distribution  in  an  infinitely  large  elastic  plate  with  a 
traction-free  elliptical  hole  under  the  state  of  simple 
tension  at  infinity  in  the  direction  of  the  X-axis  is 
presented. 

2 . 1 Formulation  of  the  Problem 
In  this  section  all  the  fundamental  elasticity  equa- 
tions to  formulate  the  so-called  "first  fundamental  problem" 
for  anisotropic  generalized  plane  stress  will  be  stated. 

We  start  from  the  generalized  elasticity  law  (Hooke's 
law)  for  orthotropic  media,  referred  to  the  symmetry  axes. 
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with  four  independent  elastic  constants  S , S , S , and 

11  E 22 

S , as  mentioned  before,  where 
66 


S = 1/E  , S = -V  /E  . 
11  1 12  12  1 


S = 1/E  , S = 1/G 
22  2 66  12 


For  the^case  that  the  fiber  direction  does  not  -coincide 
with  the  geometric  axes  as  shown  in  Figure  2,1,  the  follow- 
ing transformation  is  considered: 


Figure  2.1:  Fiber  direction  is  inclined 

to  x-axis  by  an  angle  0 

The  transformation  for  the  stresses  between  x-y  axes  and  the 
1-2  axes  will  be 
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X 

a 

y 

T 

xyj 

COS^  0 

sin^  0 

-2  sin  0 

cos  0 

sin^  0 

COS^  0 

2 sin  0 

cos  0 

sin  0 cos  0 

-sin0  cos  0 

COS^  0 - 

sin^  0 

a ■ (2.2) 

2 


T 

12 


Similarly , the  transformation  for  the  strains  .is  given  by 
the  expression 


cos^0  sin^  .0  -2sin.O  cos  0 


e 

X 

" e 

y 

1 

2^xv 

\ ^ y 

sin  0 cos^  O 2sin0  cos  0 

sin0  cosB  -sin9cos0  cos^0-sin^0 


1 

^ 12 


(2.3) 


Substitution  of  Equations  (2.2)  and  (2.3)  into  Equa- 
tion (2.1)  gives  the  generalized  plane-stress  Hooke's  Law 
referred  to  arbitrary  X-  and  Y-axes, 
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(2.4) 


Notice  that  the  coefficients  of  the  stress-strain 
matrix  in  this  case  are  all  nonzero,  which  is  similar  to 
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a generally  anisotropic  case,  but  these  coefficients  are  not 
independent.  The  expressions  of  these  coefficients  in  terms 
of  the  four  basic  constants  will  be  given  in  Appendix  A. 

The  stress  components  in  a continuous  elastic  mediiam 
in  a state  of  equilibrium,  with  no  body  force,  satisfy  the 
equations  of  equilibrium. 


9a  9t 

EL  + EiZ 

9x  9y 


0 


(2.5) 


In  the  case  of  small  strains,  the  strain-displacement 
relations  are 
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e 
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9u 

y , 

9y 


e 

xy 


9u 

X 

9y 


9u 


9x 


(2.6) 


With  Equations  (2.6)  the  displacement  u^  and  u^  can  be  elim- 
inated to  obtain  the  compatibility  equation  in  terms  of  e , 

e , and  y 
y xy 


9^e  9^e  9^e 

^ ^ - 2 ^ = 0 

3^y  9^x  9y  3x 


(2.7) 
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Substitution  of  Equation  (2.4)  into  Equation  (2,7)  gives 


g2  f 

ao+a  a+a  t 

gy2  ( 11  X 12  y 16  xy 


+ 


2 


8 

8x^ 


o + a a+a  t 
X 22  y 26  xyj 


8^ 

8x  8y 


a 

16 


o + a a+a  t 
X 26  y 66  xyJ 


0 (2.8) 


According  to  the  theory  of  elasticity,  any  state  of 
stress  must  satisfy  equilibrium  equations,  compatibility 
equation,  and  boundary  conditions.  So,  by  making  use  of  the 
Airy^  stress  function  F(x,y), 


a 

X 


8^F 
8x  8y 


(2.9) 


which  automatically  satisfies  the  equilibrium  equations, 
and  then  substituting  F(x,y)  into  the  compatibility  equation 
(2.8),  we  obtain 


a 

22 


8x“ 


2a 


8^F 


26 


8x^  8y' 


2a 


+ a 


12 


66 


8‘*F 

8x^  8y^ 


- 2a 


8'*F 


“ 8x^  8y2 


+a  i^=0 
“ 8y“ 


(2.10) 


^As  first  noted  by  a British  astronomer  G.  B.  Airy  in  1862. 
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which  can  be  written  in  the  following  sjmibolic  form 


D D D D F(x,v)  = 0 

1 2 3 <4  ' 


where 


k = 1.2,3.“ 


and  y are  the  roots  of  the  characteristic  equation 


a y^  - 2a  y^+ 
11  16 


2a  + a 
12  66 


y^-2a  y+a  =0 
26  22 


(2.11) 


On  the  basis  of  energy  considerations,  Lekhnitskii 
[18]  showed  that  the  roots  of  the  above  equation  fall  into 
one  of  the  following  cases : 

Case  1 

The  roots  are  pairwise  different 

y = a*  + i6*  and  y = a*  + iB*  (2.12a) 

111  222 


Case  2 

The  roots  are  pairwise  equal 

y = y = y = a*  + iB*  (2.12b) 

12 

where  a*,  a*.  B*.  and  B*  are  real  constants  and  i = 

12  1 2 
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In  general,  the  stress  function  F(x,y)  can  be 

expressed  by  two  arbitrary  analytical  functions  (z  ) , 

4)  (z  ) as  follows: 

2 2 

Case  1 

We  have  z =x+yy,  z =x+yy 
1 12  2 


F(x,y)  = 2 Re[(f.  (z  ) + (^  (z  )]  (2.13a) 

11  2 2 

Case  2 

We  have  z = z = x + yy 
1 2 

F(x,y)  = 2 Re[4)  (z  ) + Z(j)  (z  )]  (2.13b) 

11  2 2 


where  Re  denotes  the  real  part  of  the  complex  variables  and 
the  superbar  ( ) indicates  complex  conjugate.  Now  let 


3z, 


k=l,2 


(2.13c) 


Substituting  (2.13a)  into  (2.9),  we  get 


CT^(x.y)  = 2Re 


3$  (z  ) 9$  (z  ) 


1 1 


9z 


+ — Iz ' 


Oy(x,y)  = 2Re 


'3<I>  (z  ) 9$  (z  )■ 


1 1 


2 2 


9z 


9z 


(2.14) 


T (x,y)  = - 2Re 
xy  ' 


'9$  (z  ) 9$  (z  ) 
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Substitution  of  Equation  (2.14)  into  (2.4)  and  integration 
of  the  resulting  equation  leads  to  the  following  expression 
for  the  displacements  u and  v in  x and  y directions , 
respectively, 

u = 2Re[P  $(z)+P$(z)]-y3  +6 

111  2 2 2 ^ r u 

(2.15) 

V = 2Re[r  $(z)+r$(z)]+xB  +B 

111  2 2 2 r V 

where  the  quantities  P^ , r^  (j=*,2)  are  defined  by 


=a  y?+a  -a  u. 

1 1 J 12  1 6 3 

3 = 1 , 2 

r . 
3 

= (a  y?  + a - a y.)/y. 

12  3 22  26  3 3 

CM 

II 

and  the  additional  terms  (-yB^  + B^^)  and  (xB^  + B^)  of 
Equation  (2.15)  are  expressions  of  the  rigid  body  motion  of 
the  body.  The  coefficient,  B^ , is  linked  to  the  rotation 
of  the  body,  while  and  B^  are  linked  to  the  translation  of 
the  body  in  x and  y directions,  respectively. 

2.2  The  Boundary  Conditions 
The  stress  function  F(x,y)  must  satisfy  certain 
boundary  conditions  on  the  contour  of  the  investigated  area. 
Let  the  component  of  the  traction  vector  in  the  x-direction 
be  and  in  the  y- direct  ion  be  as  shown  in  Figure  2.2,  then 
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Figure  2.2:  Multiply  connected  region 


X = a cos  6+  T sin  6 
n X xy 


Y = T cos  6+0  sin  0 
n xy  y 


From  Figure  2.2,  we  get 

cos  6 = dy/dS  , 


sin'0=  -dx/dS 


(2.17) 


Substituting  the  above  values  for  sin  6 and  dos  6 into  Equation 
(2.17)  and  using  Equation  (2.9),  we  get 


= 9 dy  _j_  9 ^F  dx  _ d 


3F 


Y = - 


n 


3y^  dS 

3x  3y  dS 

dS  3y 

3^F 

dy  3^F  dx 

d 

. * 
3F 

3x  3y 

dS  3x^  dS 

dS 

3x 

1 J 

(2.18) 


23 


By  integrating  Equations  (2.18)  with  respect  to  S,  we  get 


9x 


s 

* 

Y dS  + c , 
n 1 

0 


3y 


X dS  + c 
Jo  n 2 

(2.19) 


where  and  are  arbitrary  real  constants.  By  substituting 
Equation  (2.13)  into  Equation  (2.18),  we  get 


2Re[$  (z  ) + $ (z  )]  = + 
11  2 2 


2Re[y  <I>  (z  ) + y $ (z  ) ] 
111  222 


s 

* 

Y dS  + c 
n n 1 


f 

1 


(2.20) 


dS  + c 

2 


f 

2 


where  the  upper  sign  corresponds  to  the  outer  boundary  L and 

2 

the  lower  sign  corresponds  to  the  inner  boundary  L , and  the 
functions  are  as  defined  in  Equation  (2.13c). 


2.3  Solution  of  the  Problem 

In  this  section  we  shall  discuss  the  solution  of  the 

problem  of  an  elliptical  hole  in  an  anisotropic  plate  under 

unidirectional  tensile  stress.  Ue  have  seen  from  Equation 

(2.20)  that  the  solution  of  the  fundamental  equation  (2.10) 

of  the  theory  of  elasticity  in  a region  which  given  traction 

boundary  conditions  on  the  boundary  of  the  region  was  reduced 

to  the  determination  of  two  analytical  functions  $ (z  ) and 

1 1 

^ (z  ) in  the  region  for  given  boundary  values  of  $ and  $ 

^ 12 

on  the  boundary  of  the  investigated  area. 
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Let  us  assume  the  two  analytical  functions  to  be  in 

the  form 


$(z)=Az  + d)(z) 

11  11  11 

(2.21) 

$(z)=Az  + d)(z) 

2 2 2 2 2 2 

where  A = E*  + iC*  and  A = B*  + iC*  are  two  complex  con- 
1 1 1 2 2 2 ^ 

stants  to  be  determined  from  the  boundary  condition  at 

O O 

infinity;  and  4>  (z  ) , 4>  (z  ) are  two  single-valued  functions 

112  2 

which  can  be  expressed  in  the  form  of  Laurent's  series. 

We  insert  Equations  (2.21)  into  Equations  (2.14)  and 
evaluate  at  infinity,  where  the  components  of  the  external 
stresses  are 


T 

xy 


0 


(2.22) 


By  carrying  out  this  substitution,  we  obtain  the  following 
linear  system  of  equations  (three  equations  with  four 
unknowns) : 


2ReIy2(B*  + iC*)  + + ic*) ] = p 

11  1 2 2 2 

2Re[(B*  + iC*)  + (B*  + iC*)]  = 0 (2.23) 

112  2 

-2Re[y  (B*  + iC*)  + y (B*  + iC*) ] = 0 
11  1 2 2 2 

In  order  to  find  a solution,  one  of  the  unknown  B*  C*  B* 

1 ’ 1 ’ 2 ’ 

and  has  to  be  chosen  arbitrarily.  So,  let  us  choose 
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C*  = 0 and  solve 
1 


B*  = 


1 


B* 

2 


c* 

2 


the  system  of  Equations  (2,23)  to  get 

P/{2[(a*  - a*)2  + (6*2  - 6*")]} 

2 1 2 2 

-P/{2[(a*  - a*)  + (6*^  - 6*")]}  = - B* 

2 1 2 1 1 

(2.24) 


P(a*  - a*)/{26*[(a*  - a*) ^ 
1 2 2 2 1 


+ (6*2  - 6*")]}  = (a*  - a*)B*/6* 

2 1 12  12 

Thus,  Equation  (2.21)  becomes 

$ (z  ) = B*z  + (^“(z  ) 

11  11  11 

(2.25) 

$ (z  ) = B*w  z + d)  (z  ) 

2 2 1 2 2 2 

where 

a*  - a* 
w = - 1 + i — 1 2- 

6* 

2 

After  substituting  Equation  (2.25)  into  (2.20), 
evaluating  at  the  inner  boundary,  and  taking  into  considera- 
tion that  the  contour  of  the  hole  is  stress  free,  i.e., 

f = f = 0,  we  can  write 
1 2 
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2Re[(|)‘*(z  ) + (})°(2  )]  = f^Ce) 

11  2 2 1 

(2.26) 

2Re[y  <j)  (z  ) + y (((“(z  )]  = f**(e) 

111  2 2 2 2 


where 


f (e)  = - B*2Re[z  + w z ] 

1 112 

(2.27) 

f (0)  = - B*2Re[y  z + y w z ] 

2 1 112  2 


After  substituting  Equations  (2.25)  into  the  stress 

formulas  of  Equations  (2.14),  using  the  value  of  B*  from 

1 

Equation  (2.24)  and  simplifying,  we  obtain  the  following 

formulas  for  a , a . and  t 

X y’  xy 


a 

X 


= P + 2Re 


94)  (z  ) 
1 1 

9z 


+ y 


(z,) 

2 2 2 


9z 


o 

y 


2Re 


94>j  (Zj) 
9z 

1 


+ 


94>  (z  ) 
2 2 

9z 

2 


(2.28) 


T 

xy 


2Re 


94>j  (z^) 
9z 

1 


+ y. 


94>  (z  ) 
2 2 

9z 

2 


By  making  use  of  the  well-known  Schwartz  formula 
[19],  which  determines  a holomorphic  function  inside  a 
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circular  region  in  terms  of  a continuous  real  function  on  its 
boundary,  then  a solution  of  Equation  (2.25)  can  be  deter- 
mined. Using  Muskhelishvili ' s notation  [19],  we  can  write 
the  Schwartz  formula  as 


F(C) 


U(e) 


g + C da 
o - C g 


(2.29) 


where  U(0)  is  the  real  part  of  F(c)  evaluated  on  the  contour 

Y of  the  unit  circle  |c|  = 1,  and  a is  some  real  constant. 

0 

To  use  this  formula  in  the  case  of  the  elliptical 
hole,  a mapping  function  has  to  be  introduced.  It  is  known 
from  the  theory  of  functions  of  a complex  variable  that  the 
mapping  of  an  area  A,  consisting  of  the  outside  of  an  ellip- 
tical hole,  onto  the  inside  of  a unit  circle,  can  be  obtained 
by  means  of  the  function 


z = o)(c) 


- b)c  + (a  + b)| 


(2.30) 


where  z-x  + iy,  C = 5 + iri,  and  a,b  are  the  semi-axes  of 
the  ellipse. 

We  introduce  two  additional  complex  planes  z and  z 

1 2 

related  to  the  z-plane  by  affine  transformation  z = x 

1 

+V^y,  z^  =x+v^y,  where  and  y^  are  as  defined  in  Equa- 
tion (2.12a).  The  known  mapping  function  z = u)(c)  then  formu- 


lates the  functions  which  give  a conformal  mapping  of  the 
z^-  and  z^ -planes  onto  the  interior  of  the  unit  circle. 
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2-Plane 


2 -Plane 
1 


2 -Plane 
2 


c;-Plane 


C = C 


The  planes  2j  and  are  obtained  from  2-plane 
by  means  of  an  affine  transformation  where 


for  the  2 -plane:  x = x + a*y, 

1 1 i'' 

for  the  2 -plane;  x = x + a*y, 

2 2 2-' 


= 3*y 
y^  = ejy 


y = a*  + iB*  y = a*  + i$* 

11  1 2 2-2 


Figure  2.3:  Mapping  the  outside  of  an  elliptical  hole, 

on  the  inside  of  a unit  circle 
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Z = (0  (c)  = 
1 1 


1 


(a  + iy  b)c  + (a  - iy  b)-^ 
1 1 C 


(2.31) 


(a  + iy  b)c  + (a  - iy  b)^ 

2 2 C 

Substituting  Equation  (2.31)  into  Equation  (2,27),  we 
obtain  the  expressions  for  the  continuous  real  functions 
evaluated  on  the  contour  of  the  unit  circle 


z — w(?)— y 
2 2 ^ 


where 


f (e)  = - B*2Re[Z 


1 


1 


+ w Z ] 
2 


f (6)  = - B*2Re[y 
2 1 


z 

1 


+ wy  Z ] 
2 2 


Z = i(a  + iy  b)a  + i(a  - iy  b)- 

1 ^ 1 ^ 1 o 

Z = i(a  + iy  b)a  + i(a  - iy  b)- 

2 ^ 2 ^ 2 o 


(2.32) 


Using  the  complex  conjugate  method  to  separate  the 
real  and  imaginary  parts  of  Equation  (2.31)  and  carrying 
out  algebraic  details,  we  obtain  the  following  expression: 


f (e)  = 
1 


4 B*k  /a 
^ 1 1 


(2.33) 


f (6)  = 
2 


~ ^ bcri  + k /a 


where  k and  k are  two.  complex  constants. 
12 
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By  using  the  Schwartz  formula,  Equation  (2,29),  and 
taking  into  consideration  that 


a 


y 


g + ^ dg 
g - C g 


A-rriC  ; 


[ — o C dg 
j o a - i g 
Y 


0 


we  obtain,  with  the  help  of  Equations  (2.13),  (2.26),  and 

(2.29),  expressions  for  the  two  analytical  functions  (|)'*,(j>° 

1 2 

in  the  ^ -plane: 


' - 2(p  u ) « (2-34) 

1 2 

' (2.35) 

1 2 

Let  us  go  back  to  the  original  z^-  and  z^ -planes.  We  invert 
the  transformation  of  Equation  (2.31)  to  give 


a - iy  b 

1 

z^  + - (a"  + y^b^) 


C 


a - iy  b 

2 

^2  + 7^2  - (a"  + y^b") 


(2.36) 


(2.37) 


and  then,  inserting  for  ^ its  value  from  Equation  (2.36)  into 
Equation  (2.34),  and  the  value  from  Equation  (2.37)  into 
Equation  (2.35),  gives 
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i P b 


a - iy^b 


2(y. 


^2^  2^  + - (a^  + y^b^) 


(2.38) 


2 2 


i P b 


a - iy  b 
2 


2(y 


1 ^2^  z + /z^  - (a^  + y^b^) 

2 V 2 2 


By  differentiating  the  fvinctions  of  Equation  (2.38),  we 
obtain 


3((i  (z  ) 
1 1 

3z 

1 


i P b F* 

1 

2(y  - y ) (a  + iy  b) 

1 2 1 


3<{)  (z  ) 
2 2 

3z 

2 


i P b F* 

2 

2(y  - y ) (a  + iy  b) 

1 2 1 


(2.39) 


where 


z 

F*  = 1 ^ 

yz^  - (a^  + y^b^) 
z 

F*  = 1 — ^ . ■ 

/z^  - (a“  + y=b^) 

Substituting  Equation  (2,39)  into  the  stress  formula  (2.28), 
yields  the  general  expression'  for  the  stress  components 
in  Cartesian  coordinates 
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(x.y) 


» 

1 

y2p*  y2p* 

- 1 + b • Im 

1 

11  2 2 

y - y 

a + iyb  a + i b 

> 

_ 1 2 

1 J 

QyCx.y)  = P B • Im 


1 

f F* 

1 

F*  1 
2 

y - y 
_ 1 2 

a + iy  b 
1 

a + iy  b 
2 J 

(2.40) 


T (x,y)  = - P b 
xy  *■’ 


1 

r y F* 
1 1 

y F* 
2 2 

^ 

a + iy^b 

a + iy  b 
2 J 

where  Im  denotes  the  imaginary  part  of  a complex  variable. 
Since  it  is  assiamed  that  no  external  forces  act  on  the  con- 
tour of  the  elliptical  hole,  then  the  normal  stress  and  the 
shear  stress  on  that  contour  should  vanish.  The  only 
remaining  stress  component  in  the  normal- tangential  coordi- 
nate directions  will  be  the  internal  normal  stress  acting  on 
planes  perpendicular  to  the  tangential  direction.  To  obtain 
a clear  picture  of  the  changes  of  this  stress  due  to  aniso- 
tropy, it  is  necessary  to  transfer  the  components  of  the 
stresses  in  Equation  (2.40)  to  elliptical  coordinates  by 
using  the  following  transformation  [21] : 

aj.  = i(a  + o)+i(a  - a)  cos  2a  + t sin  2a 
^ ^ X y 2 X y xy 

= 2 (o^  °y^  ~ 7 ^°x  ■ ^y^  ■ "^xy  (2.41) 

■ I ^®x  “ "^xy 
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where  5 and  ^ are  the  normal  and  the  tangential  axes  to 
elliptic  boundary  and  e is  the  angle  between  the  x-axis  and 
the  normal  to  elliptic  boundary  as  indicated  in  Figure  2.4, 
where 


tan  a = — tan  0 
a 


(2.42) 


Figure  2.4:  Coordinate  system  for  elliptic  hole. 


To  achieve  this  objective  we  introduce  the  parametric  form 
of  the  elliptical  contour 


X = a cos  0 , y = b sin  0 (2.43) 

and  calculate  the  components  of  stresses  o . o and  t as’ 

x’  y’  xy 

functions  of  0 on  the  contour  before  using  the  transformation 
of  Equation  (2.41). 


CHAPTER  3 

ENERGY  PRINCIPLES  AND  FINITE -ELEMENT  MODEL 

The  main  purpose  of  this  chapter  is  to  review  the 
fiandamental  bases  of  many  existing  finite  element  procedures. 
It  is  worthwhile  to  examine  how  the  energy  approaches  in 
variational  principles  play  such  important  roles  in  obtain- 
ing  approximate  elasticity  solutions  to  boundary  value 
problems.  These  variational  principles  are  based  on  the 
concept  of  a certain  integral  expression  (tt)  . When  this 
expression  has  a stationary  value  (67t  = 0)  , the  remaining 
elasticity  requirements  are  satisfied,  thus  giving  a complete 
solution , 

A traditional  application  of  the  energy  principles 
to  structural  analysis  is  the  Rayleigh-Ritz  method,  in  which 
a single  trial  function  is  constructed  to  approximate  either 
the  displacements  or  stresses,  depending  upon  the  functional 
to  be  used.  Since  this  trial  function  must  span  the  entire 
structure,  it  often  becomes  lengthy  and  complicated — an 
inconvenient  aspect.  In  an  attempt  to  overcome  this  draw- 
back, the  finite  element  method  uses  several  trial  functions, 
each  usually  simple,  and  each  defined  over  only  a portion 
of  the  structure.  This  piecewise  estimation  is  the  essence 
of  the  finite  element  method.  If  the  estimation  pieces 
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finite  elements"  are  then  properly  joined,  a satisfactory 
approximation  can  be  expected. 

There  are  two  classical  basic  variational  principles 
in  structural  mechanics.  They  are  based  upon  the  following 
principles : 

1.  Potential  energy  principle 

2.  Complementary  energy  principle 

The  other  energy  principles  have  been  generated 
from  the  previous  ones,  in  the  form  of  modified  and/or  mixed 
energy  formulations.  Some  of  the  other  energy  principles 
will  be  discussed  later  in  this  chapter. 

3.1  Potential  Energy  Principle 

To  introduce  this  principle  we  start  from  the  poten- 
tial energy  functional  (Tip)  as  given  by 

= U + V (3.1) 


where  U is  the  internal  strain  energy  and  V is  the  potential 

energy  functional  of  external  work  done  by  the  applied  loads. 

Then,  a definition  of  the  principle  of  potential  energy  can 

be  given,  as  stated  by  Desai  and  Abel  [22]: 

Of  all  possible  displacement  configurations  a body 
can  assize  which  satisfy  compatibility  and  the 
constraints  or  kinematic  boundary  conditions,  the 
configuration  satisfying  equilibrium  makes  the 
potential  energy  assume  a minimum  value.  (p . 58) 
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Here  it  is  important  to  note  that  variations  of 
displacements  are  taken  while  forces  and  stresses  are 
assiamed  constant.  Thus,  671^  = 6U  + 6V  = 0,  and  it  is 
minimal  if  > 0.  In  order  to  base  a solution  on  this 

principle,  the  potential  energy  functional  is  defined  as  a 
function  of  the  displacements  and  strains  in  the  following 
fashion  (with  no  body  force) : 


^k)l 


dV 


u. 

1 


dS 


(3.2) 


where  the  volume  integral  represents  the  work  of  internal 
strain  energy  U,  and  the  variation  of  the  surface  integral 
represents  the  external  work  done  by  the  specified  surface 
traction  (T^)  (therefore  potential  energy  lost)  during  dis- 
placements 6u^ . Then,  provided  that  the  compatibility  con- 
ditions are  enforced  by  choosing  a continuous  function  for  the 
displacement  field,  the  equilibrium  conditions  are  an 
automatic  consequence  of  the  stationary  condition  617^  = 0. 

The  finite  element  model  formed  from  this  theorem  is 
identified  as  the  conventional  finite  element  displacement 
model.  It  is  considered  the  most  predominant  finite  element 
model  to  have  been  used  in  structure  analysis.  It  provides 
a simple  yet  accurate  method  for  many  engineering  applications. 

This  model  will  be  used  and  referred  to  as  the  regu- 
lar finite  element  model  when  we  apply  it  to  the  fracture 
mechanics  aspect  of  a general  elasticity  problem  in 
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Chapter  4.  The  niathematical  derivation  of  that  regular 
finite  element  formulation  will  be  given  in  Appendix  B. 

3.2  Complementary  Energy  Principle 

Analogous  to  the  potential  energy  theorem,  the  total 

complementary  energy  functional  (tt  ) is  defined  as 

c 

= U*  + V*  (3.3) 

where  U*  is  the  complementary  strain  energy  and  V*  is  the 

complementary  energy  functional  of  boundary  forces  acting 

through  prescribed  displacements.  The  principle  can  be 

stated  as  in  Desai  and  Abel  [22]: 

Of  all  possible  force  and  stress  states  which 
satisfy  equilibrium  and  the  stress  bomdary 
conditions,  the  state  satisfying  compatibility 
makes  the  complementary  energy  assume  a minimum 
value.  (p.  59) 

Thus,  6tt  = 6U*  + 6V*  = 0,  and  ir  is  a minimum  if  >0. 

In  a manner  similar  to  that  of  potential  energy,  the 
complementary  energy  function  is  defined  as  a function  of 
stresses  and  tractions  while  displacements  are  assumed 
constant. 

The  total  complementary  energy  is  expressed  as 

"c  = 1 ? “ij  “m  ■ I ^1 

V s 

u 

In  this  expression  the  volume  integral  contains  the  comple- 
mentary strain  energy  (strain  energy  expressed  as  a function 
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of  the  stresses)  while  the  variation  of  the  second  integral 
represents  the  external  complementary  work  done  by  the 
tractions  on  the  specified  displacements.  As  long  as  the 
equilibrium  conditions  are  ensured  (a_  j = 0 for  no  body 
forces)  the  compatibility  conditions  are  implicitly 
contained  in  the  statements  of  stationary  condition  67t^  = 0. 
The  finite  element  model  based  on  this  theorem  has  been 
classified  and  identified  as  an  equilibrium  model. 

Here  the  primary  variable  is  the  stress  field,  which 
is  assumed  to  satisfy  exactly  the  equilibrium  conditions, 
while  compatibility  is  approximated  by  the  variational 
principle;  this  leads  to  better  results  for  stresses  than 
for  displacements.  Since  for  many  engineering  purposes  a 
knowledge  of  stress  distribution  is  the  desired  result  and  the 
displacements  are  of  secondary  interest,  this  model  has 
advantages  in  engineering  applications.  However,  this  com- 
plementary energy  method  has  had  limited  success  compared 
with  the  potential  energy  method.  This  limited  success  is 
due  to  the  difficulty  in  establishing  an  assumed  stress 
distribution  which  exactly  satisfies  the  equilibrium  condi- 
tions and  which  is  still  expressible  in  terms  of  a convenient 
set  of  nodal  quantities, 

3.3  Modified  Energy  Principles 
and  Hybrid  Models 

The  finite  element  procedures  based  upon  conventional 
energy  principles  (displacement  and  equilibriiam  models) 
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have  met  some  difficulties  [23]  when  elements  of  different 
types  (bars,  beams,  plates)  or  different  forms  (second  and 
third  order  polynomials)  are  joined  together.  These  diffi- 
culties have  been  encountered  in  satisfying  the  required 
conditions  on  interelement  continuity.  For  instance,  in  a 
crack  problem  the  elasticity  solution  near  the  crack  tip 
indicates  a region  of  high  stress  gradients,  with  the  stress 
becoming  singular  at  the  crack  tip.  This  known  solution 
can  be  used  to  construct  special  functions  for  elements 
around  the  crack  tip.  Meanwhile,  regular  polynomials  are 
used  for  the  other  elements  in  regions  where  the  stress 
gradients  are  not  so  high.  As  evident  in  Figure  3.1,  when 
a singular  and  a regular  element  are  adjacent,  discontinui- 
ties are  bound  to  occur  between  the  elements. 

A basic  extension  of  the  compatible  and  equilibriijm 
models  is  needed  to  handle  such  cases.  In  order  to  simply 
construct  the  different  trial  functions,  it  is  necessary 
to  relax  the  troublesome  continuity  requirements.  In 
removing  these  requirements  from  their  explicit  status, 
some  action  must  be  taken  to  reinstate  continuity  at  least 
by  some  approximate  matching  procedure.  This  can  be  ac- 
complished by  introducing  the  modified  energy  principles. 
Several  matching  procedures,  between  the  bo\mdary  of  the 
singular  cracked  and  the  adjacent  regular  elements,  have 
been  discussed  and  summarized  by  Spiering  and  de  Pater  [24] . 
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The  Lagrange  multiplier  concept  is  considered  as 
one  of  the  possibilities  to  deal  with  such  constraints  in 
variational  problem.  To  ascertain  its  usefulness,  assume 
that  a functional  (tt)  is  subject  to  the  constraint  condi- 
tion: 

on  S:  <|)(x^)  = 0 

An  augmented  functional  can  then  be  formed: 

%ug  = ^ - j \ dS  (3.5) 

s 

where  the  Lagrange  multipliers  (X^)  are  functions  of  the 
surface  coordinates.  Once  the  variation  of  the  augmented 
functional  has  been  carried  out,  the  physical  definition  of 
the  multiplier  can  be  determined. 

This  new  approach  in  finite  element  method,  which 
involves  independent  approximations  of  interior  and  boundary 
variables,  is  generally  called  the  hybrid  formulation. 

3.3.1  Modified  Complementary  Energy 
Principle  (Hybrid  Stress^ 

The  modified  complementary  energy  principle  (tt  ) 

c ' m 

can  be  formulated  by  introducing  the  traction  reciprocity 
conditions  as  equations  of  constraint  and  employing  the 
Lagrange  multiplier  concept;  that  is,  in  the  finite  element 
method,  the  equilibrium  must  be  maintained  for  the  surface 
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tractions  T , defined  by  T.  = a . , v.  across  the  inter- 
i i ij  J 

element  boundaries . Let  two  neighboring  elements  s and  r 
be  isolated  (Figure  3.1b)  and  consider  the  boxmdary  trac- 
tion components  T®(S)  and  T^(S)  (1=1,2)  over  the  respective 
sides  of  the  common  boundary  AB  for  the  two  dimension  case. 
The  equilibrium  conditions  at  the  interelement  boundary 
are  given  by 

T®(S)  + T^(S)  = 0 (1=1,2)  (3.6) 

Equation  (3.6)  may  be  considered  as  the  condition  of  the 
constraint  and  can  be  introduced  by  including  Lagrange  mul- 
tiplier terms 


X^(S) [T®(S)  + T^(S)]  dS 
AB 


(3.7) 


or 


T,  dS 

+ 

J i 

1 

CO 

AB  aB 


dS 


The  Lagrange  multipliers  X^,  which  are  functions  of  the 

surface  coordinates,  are  to  be  treated  as  additional  variables. 

When  all  the  interelement  boundaries  in  each  element 

n have  been  considered,  the  modified  complementary  energy 

functional  may  be  written  as  follows,  where  S is  the  sum  of 

n 

All  the  parts  of  the  element  boundary  which  are  interelement 
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S 

u = a + a/r  + ar 

1 2 3 


Regular  element  (r) 
u^  = b +bs+bs^ 

1 2 3 


a.  Inter-element  boundary  displacement  for 
the  assumed  displacement  hybrid  model 


B B 

T^(s)  ' 


T^(s) 

^ s 

/I 

A ■ ' 

Singular  element  (s)  Regular  element  (r) 


T®(s)  + T^(s)  = 0 i=i,2 


b.  Inter-element  boundary  traction  for 
the  assumed  stress  hybrid  model 


Figure  3.1:  Matching  regular  and  singular  elements 
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bomdaries  and  is  the  part  of  the  element  boundary  that 

n 

an  external  boundary  part  of  the  element  boundary  on  which 
displacements  u^  are  specified 


IS 


M 


" I 


n=l 


I °ij  °kn  - 


n 


dS 


n 


S T.  u.  dS 
^n  ^ ^ 


(3.8) 


Then,  by  taking  the  variation  of  the  functional  (tt  ) with 

cm 

respect  to  and  X^,  this  variation  enforces  equilibrium 
and  defines  the  multipliers  X^  as  the  displacements  u^ 
along  the  interelement  boundary.  This  modified  principle 
then  becomes  a two-field  principle,  that  is,  an  assumed 
stress  field  within  the  element  and  an  independent  assumed 
displacement  on  its  boundary. 

Finally,  since  the  hybrid  method  relies  on  assumed 
boundary  displacements,  the  prescribed  boundary  stress  no 
longer  constitutes  a restrained  boundary  condition.  In 
this  case  the  functional  can  be  written  more  con- 

veniently as  [25] 


M 


= I 

cm  ^ 


n=l 


I ‘=ijk«  °ij  “kt  ^ - 


T.  u.  dS 
1 1 


n 


8A 


n 


T . u.  dS 
1 1 


(3.9) 


n 
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where  9A^  ~ ^ entire  boundary  for  the 

element  number  n,  and 

= u.  on  (3.10) 

n 

and 


T 


i 


on  (3.11) 

n 


3.3.2  Modified  Potential  Enengv 
Principle  (Hybrid  Displacement) 

In  applying  the  variational  principle  of  minimum 
potential  energy  to  finite  element  analysis,  the  assxmed 
displacement  function  needs  to  be  continuous  within  each 
individual  element  A^,  but  not  across  the  interelement 
boundary  of  A^.  This  interelement  boundary  discontinuity 
(see  Figure  3.1a)  can  be  relaxed  by  introducing  the  Lagrange 
multipliers,  X^,  as  a constraint  condition  on  the  interele- 
ment boundary  displacement.  The  potential  energy  functional 
is  therefore  modified  and  a new  functional  can  be  written 


M 


(tt  ) = I 

P ” nil 


ifp  - I Xj  (S)I(up'  - (up'' 


] dS 


(3.12) 


AB 


Following  the  same  line  of  thought,  the  variational 


P^^^ciple  enforces  compatibility  and  defines  the  physical 
meaning  of  the  multiplier  (X^)  as  the  interelement  boundary 
tractions  T^  [26]. 
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M 


= I 


p m 


n=l 


Ti  “i  dS 


8 A 


n 


+ 


u. 

1 


dS 


(3.13) 


where  represent  the  smooth  continuous  displacements  with- 
in and  the  displacements  u^  represent  the  boundary  dis- 
placements on  the  (the  total  boundary  of  A^) . 

An  important  characteristic  of  the  displacement 
hybrid  variational  functional  is  that  the  boundary  displace- 
ment field  u^  and  the  displacement  field  u^  in  the  interior 
of  the  body  are  separate  and  distinct  [27] . 


3.3.3  Reissner's  Variational 
Principle  (Mixed  Method) 

Reissner's  functional  (tt^^)  is  one  of  the  known  mixed 
variational  principles.  It  attempts  to  treat  both  compati- 
bility and  equilibrium  conditions  with  equal  respect.  This 
functional  is  developed  by  noting  that  if  the  stresses  and 
the  strains  satisfy  the  constitutive  relations,  then 


1 


E,  . dV 
ij 


I 1 ^ijk£ 

V 


J ^ijkil  ®ij  "^ki. 


dV 


(3.14) 
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Using  this  relation,  the  complementary  energy  (in  the  form 
of  stress  energy)  can  be  brought  into  the  potential  energy 
functional : 


pc 


“ij  O'' 


■ 1 ? 


? SjtJl  “ij  “icll 


u^  dS 


(3.15) 

The  displacements  are  still  the  only  variables  since  the 
strains  and  the  stresses  are  determined  from  them.  The 
stresses  become  the  second  set  of  variables  by  relaxing  the 
compatibility  conditions  both  along  the  boundary  and  in  the 
interior  using  Lagrange  multipliers.  Introducing  these  as 
constraints,  the  functional  can  then  be  modified  [22] 


IT 


R 


ir 

pc 


- I T.  [u^  - u^]  dS  (3.16) 

s 

u 

where  the  Lagrange  multipliers,  and  T^,  become  the 
actual  stresses  and  boundary  tractions,  respectively.  Here 
the  constitutive  relations  are  the  only  equations  of  elastici- 
ty that  have  been  satisfied  exactly;  all  others  are  approxi- 
mated by  the  variational  procedure. 
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3.3.4  Mixed  Modified  Energy 
Principle  (Day's  ModeXy 


This  type  of  the  modified  energy  principles  has  been 


introduced  by  Day  [22]  in  1981.  He  employed  both  the  modi- 
fied potential  and  modified  complementary  energy  principles 
simultaneously,  so  the  displacements  and  stresses  in  the 
continuvim  can  be  treated  in  a more  equal  basis.  That  was 
accomplished  by  describing  the  continuiom  with  two  inde- 
pendent fields,  a compatible  displacement  field,  and  an 
equilibrating  stress  field.  Then  he  determined  some 
additional  terms  to  be  added  to  the  functional  in  order  to 
treat  the  problem  of  interelement  boundary  discontinuity  in 
both  displacements  and  tractions.  The  final  form  of  this 
functional  becomes 


M 


V 


V 


n 


n 


S 


S 


a 


u 


n 


n 


+ 


T,  u.  dS 
i 1 


(3.17) 


s 


n 


For  more  details  and  the  proof  of  this  method  see  refer- 
ence [23]  . 


A8 


3.4  Summary  and  Comparison 
The  various  finite  element  formulations  discussed 
above  in  this  chapter  will  be  summarized  in  Table  3.1  as 
given  by  Day  123] . 

The  main  differences  lie  in  the  manner  in  which  the 
governing  equations  of  elasticity  are  handled.  They  may 
be  either  exactly  satisfied  by  the  assumed  trial  functions 
or  may  be  approximated  by  use  of  the  variational  principle. 
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CHAPTER  4 

DEVELOPMENT  OF  A STRESS -HYBRID 
SINGULARITY  SUPER  FINITE  ELEMENT 

In  Chapter  3 we  discussed  some  details  of  the 
various  aspects  of  the  energy  principles  and  their  applica- 
tion in  forming  different  types  of  finite  element  models. 

In  this  chapter,  the  stress-hybrid  model  will  be  extended 
and  developed  in  order  to  treat  the  fracture  behavior  in 
composite  materials  and,  in  general,  anisotropic  materials. 

4.1  Brief  Outline  of  the  Method  Used 

In  the  analysis  that  follows,  the  development  of 
a stress-hybrid  singularity  super-element,  which  can  be 
utilized  to  model  the  behavior  of  a sharp  crack  in  a recti- 
linearly  (homogeneous)  anisotropic  plane,  will  be  undertaken. 

The  concept  of  super-element  was  first  developed  in 
1973  by  Tong,  Pian,  and  Larsy  [28]  for  isotropic  materials. 
Later,  in  1977,  Tong  [29]  extended  the  application  to  ortho- 
tropic material  in  which  the  axes  of  orthotropic  symmetry  are 
either  parallel  to  or  perpendicular  to  the  geometrical  axes 
of  the  plate.  Here,  in  this  study,  we  expanded  the  concept 
of  super-element  further  to  the  case  of  an  orthotropic 
material  with  its  axes  of  material  symmetry  making  some 
angle  0 with  one  of  the  geometrical  axes  of  the  plate. 
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In  Section  4.2  a full  explanation  of  the  mathe- 
matical formulation  is  given.  The  stress  functions  <j>  (^  ) 


and  of  Section  3.2  are  expressed  in  a slightly 

different  way  from  that  which  is  stated  in  reference  [29] 
Both  the  stress  intensity  factors  in  Section  4.5  and  the 
energy  release  rate  in  Section  4.5  were  computed 
using  the  finite  element  computer  code  developed. 


leads  from  the  expression  of  Equation  (3.8)  to  the  following 
variational  functional  (in  matrix  form)  see  Appendix  D,  where 
now  ir  denotes  the  complementary  energy  fxanctional  for  the 
single  element. 


Using  the  same  notation  as  Tong,  Pian,  and  Larsy 
[28] , and  implementing  both  the  divergence  theorem  and  the 
equilibrium  condition  on  the  complementary  energy  functional 


J3A 


{T}‘^{u>  dS  - i {T}^{u}  dS 

JaA 


(4.1) 


in  which 


on  3A  (4.2) 


ijkil 


in  A 


(4.3) 


where  the  vectors  {T},  {u},  and  {u}  can  be  approximated 
by 
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{T}  = 

[R]{6) 

(4.4a) 

{u}  = 

[U]{6} 

(4.4b) 

{u}  = 

[L]{q} 

(4.4c) 

in  which  [R] , [U] , and  [L]  are  matrices  defining  the 
assumed  fiinctional  relationships;  {6}  and  {q}  are  vectors 
of  vinknown  coefficients. 

A substitution  of  Equations  (4.4)  into  Equation 
(4.1)  yields 


TT  = {6}'^[G]  {q}  - ^ 


(4.5) 


in  which 


and 


[G] 


[Rl'^lL]  dS 

• 3A 


(4.6) 


[H] 


f 

1 

2 


3A 


([R]’^[U]  + [U]'^[R])  dS 


(4.7) 


By  applying  the  stationary  condition  on  Equation 
(4.5)  one  can  determine  {0}  and  then  eliminate  {$}  from 
Equation  (4.5).  Thus, 


0 = If  = [G]{q>  - [H]{B} 
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gives 


{6}  = [H]-^  [G]{q}  (4.8) 

and  by  substituting  Equation  (4.8)  into  Equation  (4.5), 
one  obtains 

TT  = {q}'^[G]'^[H]"‘ [G]  {q}  - i {q  [ G]  ^[H]  " ‘ [H]  [H]  ’ ' [ G]{q} 
or 

TT  = J {q}'^[K]{q}  (4.9) 

where 


[K]  = [G]^[H]’^[G]  (4.10) 

is  the  element  stiffness  matrix  for  the  super-element  as 
given  by  Tong,  Pian,  and  Larsy  [28]. 

Now,  in  order  to  calculate  the  components  of  this 
super-element  stiffness  matrix  [K] , the  matrices  [G]  and  [H] 
have  to  be  obtained  first.  A special  mathematical  treatment 
based  on  the  theory  of  anisotropic  plates  in  conjunction 
with  the  theory  of  elastic  fracture  mechanics  has  to  be 
undertaken . 

In  the  next  section  this  mathematical  formulation 
presented  in  detail,  and  then  it  will  be  programmed 
later  in  the  form  of  a special  subroutine  to  compute  the 
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super-element  stiffness  matrix  [K] . A proper  connection 
of  this  subroutine  to  the  main  program  is  required  in  order 
to  produce  the  global  stiffness  matrix.  This  stiffness 

then  can  be  applied  to  the  total  assembled  structure. 

4.2  Mathematical  Formulation  of 
Super-Element  Stiffness  Matrix 

In  the  following  derivation,  the  value  of  the  complex 
variables  and  Zj  will  be  based  on  Case  1 of  Equation 
(2.12a),  Chapter  2,  where  the  roots  of  the  characteristic 
equation  (2.11)  are  pairwise  different 

z = X + i y 

2,  = X + y^y  (4. 11) 

z = X + y y 
2 2 

Since  the  nature  of  the  crack  tip  singularity  is  that 
the  displacement  field  is  proportional  to  the  square  root  of 
the  distance  measured  from  the  crack  tip,  and  the  stress 
field  is  proportional  to  the  inverse  of  this  square  root 
distance,  then  the  proper  choiceof  a mapping  function  to 
account  for  this  singularity  in  the  ^-plane  (see  Figure  4.1) 
is 

2j  = w(^j)  = (4.12) 

Now  Lekhnitskii' s representation  of  the  stress  function  [18] 
can  be  utilized  to  express  the  components  of  stresses,  dis- 
placements, and  boundary  conditions  of  Equations  (2.14), 
(2.15),  and  (2.20)  in  terms  of  thenew  transformed  ^-plane 


as  follows: 


55 


z-Plane 


a.  Sharp  crack  is  assumed  to  lie  on  the 
negative  portion  of  the  complex 
z-plane  with  crack  tip  at  z = 0 


C -Plane 


b.  The  crack  lies  on  the  imaginary  axis 
of  the  transformed  ^-plane 


Figure  4.1:  Mapping  of  the  crack  region 
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(1)  For  stresses 


a = Re[y2<},'  (5  )/^  + y2^'  )/^  ] 

X 1 1 1 1 2 2 2 2 

“y  = Re[*;(5^)/5_  + 

\y  ' - Re[v_*;(£_)/£_  + 

(2)  For  displacements  (with  no  rigid  body  motion) 

u = 2 Re  [p  ((.(?)+  p (J)  (Ol 
11  1 2 2 2 

V = 2 Re[q  ({>  (5  ) + q (f>  (C  ) ] 

111  222 


(A. 13) 


(4.14 


(3)  For  boundary  conditions 


2 Re[(J.  (C  ) + (|.  (^  )]  = f* 

11  2 2 1 

2 Rety^cJ),  (?  ) + y <t>  (C  ) ] = f* 

111  2 2 2 2 


(4.15) 


As  we  can  see  from  Equation  (4.11)  and  the  assumption 

of  the  crack  which  is  to  be  laid  on  the  negative  portion  of 

the  real  axis  of  the  z-plane  (see  Figure  4.1),  z = z = z 

1 2 

on  the  crack  surface;  therefore,  the  complex  variables  C, 

will  also  coincide  on  the  crack  surface;  otherwise  5 

and  are  found  to  be  distinct  and  can  be  calculated. 

Assuming  that  the  functions  4>  (C  ) and  (j>  (^  ) of 

11  2 2 

Equations  (4.13),  (4.14),  and  (4.15)  are  properly  analytic. 
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a simple  polynomial  can  be  used  in  approximating  their  func- 
tional representations.  The  values  of  and  4>  can  there- 

1 2 

fore  be  approximated  by 


4.,  (5,)  ' I bj  : 4/5;)  = ! b,  (4.16) 


1 1 


j=l 


j=l 


J 2 


in  which  b and  b^  are  the  complex  constants 


(nonsymmetric  case) 


/\  yv 


b.  = B.  + i6.^ 

J J J+n 


or 


(4.17) 


✓V  /N 


bj  = Bj  ; ~ (symmetric  case) 


It  will  be  assxmied  that  the  surface  of  the  crack  is 

stress  free.  Thus,  T = T =0,  and  it  follows  from  Equation 

X y 

(2.20)  that  f = f =0. 

1 2 

Substituting  Equation  (4.16)  into  Equation  (4.15) 
and  using  Equation  (A. 17)  yields 


n . . 

y Re  [b  . + b . ] = 0 

j=l  J ^ J 2 

I Re[p  b + y b 5^  ] = 0 

Ijl  2J2 


(4.18) 
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from  which  it  is  only  necessary  to  satisfy 

Re[b.e^  + b.C^]  = 0 

J 1 J 2 

(4.19) 

Re[y  b.?^  + y b = 0 
1 J 1 1 j 2 


for  each  j in  order  to  assure  satisfaction  of  the  traction 
free  condition  on  each  crack  surface. 

Returning  to  Equation  (4.11)  and  introducing  the 
polar  coordinates,  r^^,  0^^,  one  can  express  as 

Zj^  = for  k=i,2  (4.20) 

where 


6j^  = tan"^  [$Jy/(x  + a*y)] 


(4.21) 


Thus,  it  is  not  difficult  to  express  from  Equation  (4.12) 
in  terms  of  the  polar  coordinates ; hence 


^k  ^ (j6j^/2)  + i sin  (j0j^/2)] 

Along  the  crack  surface  r = |x|,  0 = -rr 

Iv 

53  = = |x|j/2[cos  (U/2)  + i sin  (jTT/2)] 

(4.22) 
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Substituting  Equation  (4.22)  into  Equation  (4.19)  gives 
Re{  (B^  + iBj^^  + 3^  + iBj^^)  |x[  j [cos  (j  7t/2)  + i sin  ( j ir/2)  ] } = 0 

Re{  [y  (B.+iB..  )+y  (B.  +iB.  .„)]  lx|  [cos  (jir/2) 

1 J 3+n  2 J j+n  ' ' 

+ i sin  (j  tt/2)  ] } = 0 
or 

(Bj  + B^)  cos  (jTT/2)  - (3j^^  + ^j+n^  (j^/2)  = 0 (4.23a) 

(a*B  - B*B  + a*B.  - B*B.  ,„)  cos  (jr/2) 

ij  ij+n  23  2j+n 


- (B*B.  + a*B..  + B*B, 

13  1 3+n  23 


+ a*B.,  ) sin  (jir/2)  = 0 

2 3+n  *’ 


(4.23b) 


By  solving  Equations  (4.23)  for  the  real  coefficient  Bj  and 

/s 

B. . in  terms  of  B.  and  B..  , we  can  find  a relation  between 
j+n  J j+n’ 

<|)^(C^)  and  guarantee  a stress-free  condition  along 

the  entire  length  of  the  crack  face.  Thus,  the  following 
relations  can  be  obtained: 
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in  which 


3.  = M B + M B 
j iij  j i2j  n+j 


B = M B + M B 

j+n  2ij  j 22j  n+j 


(4,24) 


M 


- B*/6* 
1 2 


M = 0 , 

2lj 


M 


(a*  - a*)/B* 

2 12 

for  odd  j 


M = -1 
22j 


and 


(4.25) 


M = -1 
1 1 j 


M 


M = 0 , 
12j 


for  even  j 


21  j 


- a*)/B*  . M . = - B*/6* 

1 22  223  12 


Thus,  the  two  analytical  stress  functions  of  Equation  (4.16) 
will  become 


n 


(5  ) = I (6,  + i6,.  )C^ 

11  4=1  1 j+n  1 


(4.26) 


n 

) = I KM  + M ,B  .)  + i(M  B.  + M B )]?^ 


2 2 


3 = 1 


llj  j I2j  n+j 


21 j j 22j  n+j  2 


Substituting  the  above  value  of  stress  functions  from  Equa- 
tion (4.26)  into  the  stress  equation  (4.13),  we  obtain 
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a 

X 


I 


j=l 


A,(?  )3, 

J 1 2 J 


a 

y 


211 


(4.27) 


2n 


where 


A.  = j Re[yn^"^  + y^(M  + i M 

J 11  2 1 1 j 2 1 j 2 

A. ^  =j  Reliy^j"^  + y^(M  + i M 

3+n  11  212j  22j2 

B = j Re[^j-2  + (M  + i M 

j 1 llj  21j2 

B. ^  = j Re[i^j"2  + + i M 

J+n  1 12j  22j  2 

C = - j Re[y  + y (M  + i M 

J 11  2 1 1 j 2 1 j 2 

c.^  = - j Re[iy  + y (M  + i M )^~^] 

J+n  "^11  "^2  12j  22j  2 

A similar  substitution  of  Equation  (4.26)  into  the  displace- 
ment equation  (4.14)  produces 
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2X1 


where 


u = I f. (C  )3. 
j=l  J ^ ^ J 


(4.28) 


2n 


f . (5  ,C  ) = 2 Re[p  + p (M  + i M )C^] 
212  ^11  211j  2lj2 


(C  ,5  ) = 2 Re[i  p + p (M  . + i M 
J+ni2  11  212J  2232 


g.  (C  .0  = 2 Re[q  + q (M  + i M 
J12  IX  211J  2132 


g.^  (5  .5  ) = 2 Reti  q + q (M  + i M )C^] 

3 +>^12  11  2123  22j2 


The  integrations  to  determine  [H]  and  [G]  for  the 
particular  geometry  of  Figure  4.2  become 


{3}  [HHe)  = 


{T}^{u>  dS 


9A 


C F H 


_ 

+ 

+ 

» « 

V (a  u + T v)  dS  + 

. 

. 

X X xy 

B 

D 

G 

V (t  u + a v)  dS 
y xy  y ' 


(4.29) 
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{6}  lG]{q}  = 


{T}^{u}  dS 


9A 


fC 

|-F 

fH  , 

s 

+ 

+ 

+ ■'xy^>  * 

B 

D 

G 

(4.30) 


where 


V ' K •*•  Vj>BjSk 


(=jSk  + =kSj>8j8^ 


Vljj.i  <®jSk  ®k8j>6jBk 


In  applying  Equation  (4.30)  the  boundary  displacement 
vector  {u}  will  be  expressed  in  terms  of  the  element 
nodal  displacement  vector  {q}  by  the  following  relation 


{Q}  = [L]{q} 


(4.31) 
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where  [L]  is  a matrix  of  properly  selected  shape  functions 
defined  along  the  element  boundary.  The  [L]  matrix  is 
chosen  such  that  when  the  corresponding  nodal  displacements 


Figure  4,2:  Nine-node  super-element 

of  the  adjacent  elements  are  matched,  u is  the  same  for  the 
two  different  kinds  of  elements  over  their  common  boundaries . 
In  this  case,  the  standard  elements  are  assumed  to  be  four- 
node  quadrilaterals  possessing  four  corner  nodes  for  each 
element.  Since  these  elements  possess  a linear  displacement 
along  the  element  boundaries,  it  is  necessary  to  require  a 
linear  variation  of  the  boundary  displacement  field  {u} 
[Equation  (4.31)]  along  the  sides  of  the  singular  element 
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(Figure  4.2)  in  order  to  insure  compatibility  between  the 
two  element  types.  Therefore  one  can  assume  the  following: 

^11  ~ ^22  “ ^ ~ ^24  ~ other 

L_  = 0 and  = {u^  Up_j_^  between  nodes  P and 

nodes  P+1,  where  i is  the  distance  between  the  two  nodes, 
(^pJVp)  is  the  value  of  u at  node  P,  and  S is  the  distance 
measured  from  node  P. 

4.3  Numerical  Integration 
The  evaluation  of  the  integrals  involved  in  the 
stiffness  matrices  was  accomplished  using  numerical  integra- 
tion. In  this  regard,  it  is  noted  that  the  type  and  order 
of  integration  used  represents  a critical  decision,  as  a 
method  must  be  chosen  which  uses  a minimum  number  of  inte- 
gration points  and  achieves  an  acceptable  level  of  accuracy. 
The  method  used  here  in  formulating  the  super-element  stiff- 
ness matrix  was  a 5X5  Gaussian  quadrature  procedure;  the 
method  for  the  formulation  of  the  regular  element  stiffness 
matrix  was  a 3X3  Gaussian  quadrature.  These  choices  seem 
to  offer  the  best  compromise  between  accuracy  and  economy. 
That  is,  the  stress  intensity  factor  K^.  (opening  mode)  for 
the  first  example  used  by  Tong,  Pian,  and  Larsy  [28]  was 
obtained  exactly  in  CPU  time  0.51  seconds,  compared  to 
1.7  seconds  on  an  IBM  370/165  in  that  reference.  For  the 
other  example  of  that  study  [28]  the  exact  values  of  Kj 
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and  Kjj  (shear  mode)  were  obtained  using  the  same  number 
of  degrees  of  freedom  with  only  1.15  seconds  of  CPU  time 
compared  to  2.5  seconds;  i.e. , we  established  a 100 
percent  accuracy  with  more  than  50  percent  reduction  in 
computation  time  for  both  cases  studied  in  that  reference. 

4.4  Convergence  Characteristics  of  the  Method 
In  any  acceptable  numerical  formulation,  the  numeri- 
cal solution  must  converge  and  tend  toward  the  exact  solu- 
tion of  the  problem.  For  the  finite  element  method,  the 
pure  models  (regular  stress  and  displacement  models)  are 
based  strictly  on  the  energy  theorems,  which  are  not  only 
stationary  principles  but  are  also  minimum  principles. 

This  allows  bounds  to  be  put  on  the  strain  energy.  The 
compatible  (displacement)  model  will  always  provide  an 
upper  bound  to  the  true  stiffness  of  the  structure.  In 
other  words,  the  stiffness  coefficients  for  a given  dis- 
placement model  have  magnitudes  higher  than  those  for  the 
exact  solution.  Thus,  under  a given  load,  the  simulated 
structure  deforms  less  than  the  actual  structure.  It 
follows  that  the  approximate  displacement  solution  will 
converge  to  the  exact  solution  from  below:  that  is,  we 

obtain  lower  bounds  to  the  true  solution  [22].  The  equi- 
(stress)  model  is  just  the  opposite,  always  being 
too  flexible,  and  gives  an  upper  bound  [30]. 
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In  the  hybrid  formulations, the  relaxing  of  the  continuity 
requirements,  by  introducing  the  Lagrange  multiplier,  does 
alter  the  convergence  characteristics.  As  a consequence, 
neither  monotonic  convergence  nor  determination  of  lower  or 
upper  bounds  can  be  achieved.  However,  all  hybrid  formula- 
tions possess  a solution  that  lies  between  the  limits  [30] 
such  as  shown  in  Figure  4.3. 


1.  Equilibriimi  (monotonic) 

2.  Hybrid  (nonmonotonic) 

3.  Hybrid  (nonmonotonic) 

4.  Compatible  (monotonic) 


Figure  4.3:  Convergence  of  strain  energy 

with  the  mesh  refinement 
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In  this  study  we  chose  the  interpolation  function 
[L]  of  Equation  (4.31)  to  be  of  the  linear  type  because  it 
is  well  recognized,  in  signular  elasticity  problems,  that 
the  rate  of  solution  convergence  is  generally  independent  of 
the  order  of  pol3momials  used  for  interpolation  functions  in 
element  formulation  [31].  That  is,  the  use  of  higher-order 
elements  does  not  necessarily  improve  the  rate  of  convergence 
of  the  numerical  solution.  In  fact,  the  rate  of  convergence 
of  the  finite  element  solution  is  governed  by  the  order  of 
the  singularity  rather  than  the  order  of  the  polynomial  used 
for  the  interpolation  function.  In  addition,  simplicity  and 
less  computer  time  can  be  gained. 

4.5  Stress  Intensity  Factors 
As  previously  derived  by  Sih,  Paris,  and  Irwin  [32], 
the  stress  intensity  factors,  Kj  and  for  plane  gener- 

ally anisotropic  media  near  the  crack-tip  region  can  be 
evaluated  from  the  formulas 


^I  ~ “ - y )/y  lim  ) 

1 2 2 ^ ^0 
1 


(4.32) 


and 


^l  ^ll'^'^  ~ *^(y  - y )/y  lim  4>'(c  ) 

1 12  1 ^ 2 2 


(4.33) 
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where  from  Equation  (4.9)  we  obtain 


lim  (|)'(^  ) = b and  lim  ) = b 

f -1-0  ^ ^ ^ P -»-0  2 2 1 

1 2 


and  from  Equations  (4.15a)  and  (4.15b)  when  j=l,  we  get 


B = - 3 

n+ 1 n+j 


6 = [(a*  - a*)3  - B*6  ]/B* 

1 2 1 n+1  112 


Now  Equations  (4.32)  and  (4.33)  become 


K + K /y  = - /?[(y  - y )/y  ]b 

i ii  2 12  2 1 


K + K /y  = /I[(y  - y )/y  ]b 

•*-  iJ-l  1 211 


Solving  the  above  two  equations  for  we  obtain 


K = - /^(y  b + y b ) (4.38) 

J-J-  2 111 


Substituting  Equation  (4.38)  into  Equation  (4.36),  we  get 


= /7(b  + b ) (4.39) 

ill 


Since  is  a real  quantity  (the  stress  intensity  factor  for 
the  crack  opening  mode)  only  the  real  part  of  the  right-hand 
side  of  Equation  (4.39)  will  be  considered,  i.e.. 


(4.34) 

(4.35) 

(4.36) 

(4.37) 
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K = /I  (B  + B ) (4.40 

J-  1 1 

and  by  using  Equation  (4,35)  to  eliminate  § from  Equation 
(4,40),  a direct  expression  for  can  be  obtained 

K = /?  t(l  - 3*/6*)6  + (a*  - a*)B  /6*]  (4.41) 

12  12  1 n+ 1 2 

To  evaluate  the  same  procedure  can  be  followed 

by  considering  only  the  real  part  of  Equation  (4.38)  which 
is  to  be  evaluated  after  the  replacement  of  the  roots  y and 
y^  by  their  complex  expressions  from  Equation  (2.12a),  i.e.. 

Re [(a*  + iB*) (6  + i§  ) + (a*  + iB*) (B  + iB  )] 

2 2 1 n+1  1 1 1 n+1 

= - n (ct*B  - B*B  + a*B  - B*B  ) 

21  2 n+1  11  1 n+i 

By  eliminating  B and  B and  using  Equations  (4.34)  and 

1 n+ 1 

(4.35),  a direct  expression  for  can  be  obtained 

K,  = - /7  {(a*  - a*B*/B*)B 
12  12  1 

+ [B*  - B*  + a*(a*  - a*)/B*]B  (4.42) 

2 1 22  12  n+1 

Thus , once  the  value  of  the  parameters  B and  B have  been 

1 n+1 

computed,  the  stress  intensity  factors,  in  mode  I and  mode  II, 
can  be  calculated  directly  from  the  preceding  equations . 
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In  the  case  of  isotropic  material,  the  roots  of 

the  characteristic  equation  are  equal  to  i,e., 

y = y = i;  in  that  case  it  is  clear  that  this  value  cannot 
1 2 

be  applied  directly  to  Equations  (4.32)  and  (4.33),  but  it 
must  be  approached  in  the  limit  as  a small  difference  A, 
i.e.,  when  A approaches  zero  in  which 

P = i , y = (1  + A)i 
1 2 


or 


3*  = 1 ,-  6*  = 1 + A , and  a*  = a*  = 0 (4.43) 

12  12 

IBy.. substituting  Equation  (4.43)  into  both  Equations  (4.41) 

and  (4.42)  and  taking  the  limit  on  A in  each  equation,  and 

considering  a symmetric  deformation  (3  = 0)  and  single 

n+ 1 

fracture  mode  behavior,  we  obtain 
From  Equation  (4.41):  = /2"  3^ 

From  Equation  (4.42)  : = 0 

which  are  the  well-known  formulas  of  the  stress  intensity 
factor  for  the  case  of  a single  opening  mode  in  an  isotropic 
material  [2] . 
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4.6  Strain  Energy  Release  Rate 
Although  stress  intensity  factors  and 
describe  details  of  a stress  field  near  the  crack-tip, 
the  critical  strain  energy  release  rate  at  the  begin- 

ning of  crack  propagation  is  also  of  a significant  interest. 
It  is  a quantity  physically  measurable  in  experiments  and 
methematically  well  defined.  In  addition,  it  can  be  viewed 
as  the  portion  of  the  work  performed  by  the  external  loads 
during  crack  extension  that  goes  into  failure  of  the 
material  per  unit  area  of  crack  growth. 

Sih,  Paris,  and  Irwin  [32]  studied  the  energy 
release  rates  for  the  case  of  more  than  one  mode  present  at 
crack-tip  in  generally  anisotropic  media.  They  determined 
that  this  energy  can  be  evaluated  from  the  following 
formulae : 


J = iu  K,  a Im{[KT(y  + y ) + K,,]/y  y } 

1 ^ ■1-22  -l-l  2 -1--1-  12 

3=^71  K,,  a Im[K,,(y  + y ) + K,y  y ] 

2 ^ ■L-L  1 1 -L-l-  1 2 -1-12 


(4.44) 


and 


G = J + J 
1 2 

For  the  special  case  of  orthotropic  material  with  the  crack 
being  on  one  plane  of  elastic  symmetry  Equation  (4.44)  can 
be  reduced  to  the  following  form: 
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J = ttK"  {[(a  /a 
1 2 2 11 


+ (2 


1 2 


a )/2 
6 6 


a ] 
1 1 


1 1 


a /2} 
2 2 


/2 


J 

2 


(4.45) 


ttK^  [(a  /a  + (2  a 

II  22  1 1 12 


+ a )/2  a 
6 6 11 


Where  the  values  of  a a .a  . and  a can  be  found  in 

II  22  12  66 

Appendix  A, 

4.7  Stress  and  Displacement  near  the  Crack  Tip 
After  calculating  the  stress  intensity  factors 
and  Kjj  (opening  and  shear  modes,  respectively)  from  Section 
4.5,  the  stresses  and  displacements  in  the  vicinity  of  the 
crack  tip  can  be  easily  obtained  using  the  corresponding 
mode  equations  which  are  stated  in  Appendix  C for  both 
generally  anisotropic  and  isotropic  cases. 

Here  we  notice  that  in  the  orthotropic  case  the 
stress  distributions  and  the  displacement  field  near  the 
crack  tip  depend  not  only  on  Kj  and/or  as  in  the 
isotropic  cases  but  also  depend  on  the  material  properties 
and  the  material  symmetry  axis  with  respect  to  geometric 
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CHAPTER  5 


RESULTS  OF  ANALYTICAL  AND  NUMERICAL  TECHNIQUES 

This  chapter  will  discuss  results  in  two  parts.  The 
first  will  be  the  application  of  Savin's  analytical  solution 
to  some  practical  problems  in  unidirectional  composite 
materials,  mainly  glass/epoxy  and  graphite/epoxy,  both 
widely  used  in  aircraft  and  military  industries.  The  plane 
stress  field  will  be  determined  for  thin  plates  of  these 
materials  with  the  presence  of  either  a circular  hole  or 
an  elliptical  hole.  Three  different  elliptical  holes  will 
be  discussed  in  terms  of  the  aspect  ratio  of  major  to  minor 
axis . 

The  second  part  of  this  chapter  will  discuss  the 
numerical  solution  for  the  fracture  mechanics  problem.  Stress 
intensity  factors  obtained  by  finite  element  solution  will 
be  compared  to  other  solutions  obtained  by  different 
techniques  using  the  same  assxomptions . 

5.1  Stress-Concentration  Factors 
The  problem  of  stress  concentrations  around  openings 
in  an  anisotropic  plate  was  one  of  primary  concern  to  many 
investigators  [33].  Since  the  establishment  of  the  commer- 
cial use  of  composite  materials  in  the  early  1960s,  it  has 
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become  vary  important  to  analyze  the  prediction  of  failure 
due  to  these  high  stress  concentrations. 

In  order  to  apply  any  of  the  proposed  failure  cri- 
teria for  composite  materials,  the  stress  distributions  due 
to  the  external  load  before  failure  have  to  be  determined 
first.  Because  of  the  limited  scope  of  this  study,  investi- 
gation of  these  failure  criteria  will  not  be  discussed 
here.  All  results  of  stress  distributions  around  an  opening 
to  be  given  here  are  obtained  by  applying  the  theory  of 
elasticity  according  to  the  formulation  discussed  in 
Chapter  2 under  the  assumption  that  failure  has  not  yet 
occurred. 

The  elastic  properties  used  for  calculation  of  these 
stresses  in  graphite /epoxy  are  given  here  as  in  Collazo's 
study  [34] : 

E = 128  GPa  (20  x 10®  psi) 

1 

E =11  GPa  (1.6  X 10®  psi) 

2 

G = 6.4  PGa  (0.93  x 10®  psi) 

1 2 

V =0.38 
1 2 

The  fiber  volume  fraction  is  approximately  63  percent.  The 
elastic  properties  for  the  glass/epoxy  composite  will  be 
given  in  Chapter  6,  where  this  composite  is  used  for  the 
experimental  part  of  this  study. 


76 


Stress  concentrations  will  be  discussed  for  four 
different  hole  shapes,  two  different  fiber  materials,  and 
five  different  fiber  orientation  angles  with  respect  to  the 
loading  direction.  In  all  following  figures,  which  contain 
stress-concentration  factor  vs  the  coordinate  angle  6,  we 
use  the  symbol <0  for  the  corresponding  isotropic  case,  and 
the  symbols  e,  o.  A,  +,  and  x to  stand  for  the  different 
fiber  orientation  angles  0°,  30°,  45°,  60°,  and  90°, 
respectively.  Both  the  fiber  orientation  angle  0 and  the 
coordinate  angle  0 (abscissa  in  figures)  for  calculating 
stresses  are  measured  counterclockwise  from  the  loading 
direction.  Hereafter,  we  will  discuss  the  stress-concen- 
tration factor  for  four  different  cases,  one  circular  and 
three  different  elliptical  holes. 

5.1.1  The  Circular  Hole 

The  stress  concentrations  around  the  circular  hole 
in  glass/epoxy  and  graphite /epoxy  composite  plates  subjected 
to  axial  load  (P)  are  calculated  for  different  fiber  orien- 
tation angles.  Figure  5.1  shows  the  stress  concentrations 
around  the  circular  hole  in  plates  made  of  glass/epoxy 
each  of  which  has  different  fiber  orientations,  while 
Figure  5.2  shows  the  stress  concentrations  for  the  same 
fiber  orientation  but  for  plates  made  of  graphite/epoxy 
composite.  In  both  figures  and  the  following  figures,  the 
case  of  isotropic  material  which  corresponds  to  same  hole 
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size,  shape,  and  loading  condition  is  also  presented  for 
comparison.  For  example  in  either  Figures  5.1  or  5.2,  the 
maximum  value  of  the  stress-concentration  factor  for  a 
circular  hole  in  an  isotropic  material  occurs  at  90°  and 
has  a value  of  3,  which  is  well-known  in  the  theory  of 
elasticity  [21].  Since  the  maximum  stress  concentrations 
and  their  location  angles  are  different  depending  on  the 
fiber  material  and  the  fiber  orientation  angle.  Table  5.1 
presents  a numerical  value  for  the  maximum  stress-concen- 
tration factor  (Oq/P)„„^  and  its  location  (angle  6*)  on 
the  circular  contour  of  the  hole  for  each  fiber  orientation 
angle  and  for  the  two  fiber  materials. 

The  effect  of  fiber  orientations  on  the  maximum 
stress-concentration  factor  and  on  its  location  are  shown 
in  Figure  5.3  and  5.4,  respectively,  for  glass/epoxy  and 
graphite /epoxy  composite  materials. 
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igure  5.1:  Stress-concentrations  around  a circular  hole  for  glass/epoxy 

with  different  fiber  orientations 
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Table  5.1:  Maxiimim  stress-concentration  factors  for  a circular 

hole 


Fiber 

Glass/Epoxy 

Graphite/Epoxy 

orientation 
0 (deg) 

e*  (deg) 

e (deg) 

0 

4.630 

90 

6.274 

90 

30 

3.953 

70 

5.131 

67 

45 

3.241 

63 

3.989 

57 

60 

2.626 

70 

2.870 

51 

90 

2.674 

90 

2.492 

90 

Note:  For  the  case  of  isotropic  materials  (c^ /P)^^^=  3 . 0 

and  6 =90 
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Fiber  orientation  Q (deg) 

Figure  5.3:  Maximum  stress-concentration  factor  vs 

direction  of  fiber  orientation  for  a 
circular  hole 


105 


Fiber  orientation  0 (deg) 


Figure  5.4:  Location  of  the  maximum  stress-concentration  vs 

direction  of  fiber  orientation  for  a circular 
hole 
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5.1.2  The  Elliptical  Hole  with  Aspect  Ratio  a/b  = 3 

The  stress  concentrations  around  the  hole  for  the 
two  fiber  materials  and  the  five  fiber  orientation  angles 
in  this  case  are  calculated  and  presented  in  Figure  5.5 
and  5.6  for  glass/epoxy  and  graphite /epoxy , respectively. 

Table  5.2  presents  the  n^americal  value  of  the  maximum 
stress-concentration  factor  (o_/P)_„„  and  its  location 

Tj  HlaX 

(angle  0*)  on  the  elliptical  contour  of  the  hole  for  each 
fiber  orientation  angle  and  for  the  two  fiber  materials. 

The  effect  of  fiber  orientations  on  the  maximum 
stress-concentration  factor  and  on  its  location  are  shown 
in  Figures  5.7  and  5.8,  respectively,  for  glass /epoxy  and 
graphite /epoxy  composite  materials. 


ISOTROPIC 

>e=o* 
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Figure  5.5;  Stress-concentrations  around  an  elliptical  hole,  with  ratio  a/b 
for  glass/epoxy  with  different  fiber  orientations 
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Table  5.2:  Maximum  stress-concentration  factors  for  an  elliptical 

hole  with  ratio  a/b=3.0 


Fiber 

Glass/Epoxy 

Graphite/Epoxy 

orientation 
0 (deg) 

6 (deg) 

6 (deg) 

0 

12.041 

90 

16.822 

90 

30 

9.706 

81 

13.565 

81 

45 

7.927 

79 

10.300 

76 

60 

6.014 

78 

7.016 

71 

90 

6.220 

90 

5.475 

90 

Note:  For  the  case  of  isotropic  materials  (q, /P)  =7.0 

^ M /F)3X 

and  6 = 90.0® 
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Figure  5.7;  Maximum  stress-concentration  factor  vs  direction 
of  fiber  orientation  for  an  elliptical  hole,  a/b = 3 
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0 (deg) 


Location  of  the  maximiom  stress-concentration  vs 
direction  of  fiber  orientation  for  an  elliptical 
hole,  a/b  = 3 


Figure  5.8; 
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5.1.3  The  Elliptical  Hole  with  Aspect  Ratio  a/b  = 9 

The  stress  concentrations  around  the  hole  for  the 
two  fiber  materials  and  the  five  fiber  orientation  angles 
in  this  case  are  calculated  and  presented  in  Figure  5.9 
and  5.10  for  glass/epoxy  and  graphite /epoxy , respectively. 

Table  5.3  presents  the  numerical  value  of  the  maximum 
stress-concentration  factor  (o_/P)_„„  and  its  location 

i)  UlaX 

(angle  6*)  on  the  elliptical  contour  of  the  hole  for  each 
fiber  orientation  angle  and  for  the  two  fiber  materials. 

The  effect  of  fiber  orientations  on  the  maximum 
stress-concentration  factor  and  on  its  location  are  shown 
in  Figures  5.11  and  5.12,  respectively,  for  glass/epoxy 
and  graphite /epoxy  composite  materials. 
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Figure  5.9:  Stress-concentrations  around  an  elliptical  hole,  with  ratio  a/b 

for  glass/epoxy  with  different  fiber  orientations 
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Table  5.3:  Maximiim  stress- 

hole  with  ratio 

concentration 
a/b  = 9 

factors  for 

an  elliptical 

Fiber 

Glass/Epoxy 

Graphite/Epoxy 

orientation 

* 

it 

0 (deg) 

6 (deg) 

6 (deg) 

0 

34.122 

90 

48.469 

90 

30 

27.898 

87 

39.179 

87 

45 

22.293 

86 

29.872 

85 

60 

15.632 

88 

20.184 

83 

90 

16.066 

90 

14.426 

90 

Note:  For  the  case  of  isotropic  materials  (a  /P)  = 19.0 

*0  n 

and  6 =90 
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Figure  5.11;  Maximum  stress-concentration  factor  vs  direction 
of  fiber  orientation  for  an  elliptical  hole,  a/b = 9 
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6 (deg) 


Figure  5.12;  Location  of  the  maximum  stress-concentration  vs 
direction  of  fiber  orientation  for  an  elliptical 
hole , a/b  = 9 
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5.1. A The  Elliptical  Hole  with  Aspect  Ratio  a/b  = 20 

The  stress  concentrations  around  the  hole  for  the 

two  fiber  materials  and  the  five  fiber  orientation  angles 

in  this  case  are  calculated  and  presented  in  Figures  5.13 

and  5.1A  for  glass/epoxy  and  graphite /epoxy , respectively. 

Table  5. A presents  the  mamerical  value  of  the  maximum 

stress-concentration  factor  (o  /P)  and  its  location 

n max 

(angle  6*)  on  the  elliptical  contour  of  the  hole  for  each 
fiber  orientation  angle  and  for  the  two  fiber  materials. 

The  effect  of  fiber  orientations  on  the  maximum 
stress-concentration  factor  and  on  its  location  are  shown 
in  Figure  5.15  and  5.16,  respectively,  for  glass/epoxy 
and  graphite/epoxy  composite  materials. 
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Table  5.4:  Maxinuim  stress-concentration  factor  for  an  elliptical 


hole 

with  ratio 

a/b  = 20 

Fiber 

Glass/Epoxy 

Graphite/Epoxy 

orientation 
0 (deg) 

e (deg) 

\nax 

6 (deg) 

0 

74.603 

90 

106.483 

90 

30 

59.951 

89 

72.795 

89 

45 

• 

48.576 

88 

63.407 

88 

60 

35.687 

88 

44.003 

87 

90 

34.479 

90 

30.835 

90 

Note: 


For  the  case  of  isotropic  materials 

* O 

and  6 =90 


41.0 
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Figure  5.15;  Maximum  stress-concentration  factor  vs  direction 
of  fiber  orientation  for  an  elliptical  hole,  a/b=20 
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Figure  5.16;  Location  of  the  maximiom  stress-concentration  vs 
direction  of  fiber  orientation  for  an  elliptical 
hole,  a/b  = 20 
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5.1.5  Effect  of  Aspect  Ratio  on  Stress-Concentrations 

The  maximiom  value  of  stress-concentration  factors  vs 
the  aspect  ratios  of  the  elliptical  hole  (a/b)  are  shown 
in  Figure  5.17  and  5.18  for  glass/epoxy  and  graphite /epoxy , 
respectively.  A linear  variation  of  the  maximum  stress- 
concentration  factor  and  the  aspect  ratio  has  been  noticed 
for  the  last  two  figures  for  all  fiber  orientation  angles 
as  well  as  for  the  isotropic  case.  This  linear  variation 
concludes  what  is  known  in  the  field  of  fracture  mechanics 
as  the  stress  singularity  near  crack  tip,  where  this  stress 
concentration  factor  reaches  infinity  as  the  aspect  ratio 
of  the  elliptical  hole  goes  to  infinity,  in  other  words  as 
the  minor  axes  of  the  ellipse  becomes  zero  (the  case  of 
slit  or  sharp  crack). 
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15  10  15  20 

The  ratio  a/b 

Figure  5.17:  Maximum  stress-concentration  factor  vs  aspect  ratio 

a/b  of  elliptical  hole  for  glass/epoxy 
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Figure  5.18:  Maximum  stress-concentration  factor  vs  aspect  ratio 

a/b  of  elliptical  hale  for  graphite/epoxy 
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5.2  Stress  Intensity  Factors 

After  developing  the  finite  element  computer  program 
which  is  based  on  the  stress-hybrid  model,  the  author  chose 
some  numerical  examples  which  were  based  on  different  tech- 
niques such  as  the  boundary  collocation  method.  Two 
numerical  examples  were  taken  from  Gandhi's  study  [7]  and 
analyzed  according  to  the  finite  element  procedure.  The 
first  one  has  a central  inclined  crack  with  45°  to  loading 
axis  with  the  fiber  direction  parallel  to  the  loading 
direction.  The  second  one  is  exactly  the  same  as  the  first 
one  except  that  the  fiber  direction  is  perpendicular  to  the 
loading  direction.  Figure  5 . 19  shows  a computer  mesh  of 
the  plate  with  the  singularity  elements  at  each  side  of 
the  crack  tip. 

The  stress  intensity  factors  obtained  from  the 
finite  element  computer  code  give  the  following  result: 

Kj  = 0.4561  Kjj  = 0.6223  (5.1) 

Comparing  to  the  values  obtained  from  the  boundary  collection 
method  in  reference  [7] 

Kj  = 0.6  Kjj  = 0.56  (5.2) 

we  notice  here  some  differences  which  will  be  discussed  at 
the  end  of  this  section. 

In  the  other  numerical  example  where  the  fiber 
direction  is  perpendicular  to  the  loading  direction,  the 
finite  element  result  gives 
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Figure  5.19:  Mesh  plot  before  deformation 
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Kj  = 0.3825 


K 


II 


= 0.5669 


(5.3) 


Comparing  this  to  the  values  obtained  from  reference  [7], 


indicates  also  some  difference  in  the  evaluation  of  these 
stress  intensity  factors. 

Examination  of  Gandhi's  results  [7],  Equations  (5.2) 
and  (5.4)  for  both  cases,  shows  that  the  value  of  the  stress 
intensity  factor  be  obtained  for  the  opening  mode  seems  to 
have  no  large  difference  between  these  two  cases  in  spite 
of  having  a large  difference  in  the  fiber  orientations.  For 
the  first  case  the  fiber  direction  is  parallel  to  the  load- 
ing direction  and  for  the  second  case  the  fiber  direction 
is  perpendicular  to  the  loading  direction.  In  addition, 
the  first  case  should  produce  stronger  material  in  the 
loading  direction  than  the  other  case.  That  means  stress 
intensity  factor  should  have  a higher  value  than  the 
second  case,  which  is  predicted  from  the  results  obtained 
by  our  finite  element  method.  Figure  5 . 20  indicates  a computer 
mesh  plot  of  the  plate  after  deformation.  This  figure  also 
shows  that  the  central  crack  starts  to  open  due  to  the 
external  applied  load.  These  crack  displacements  were 
magnified  about  100  times  larger  than  the  actual  displace- 
ment in  order  to  see  the  shape  of  such  deformations  on 
the  plot. 


Kj  = 0.64 


Kjj  =0.61 


(5.4) 
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Figure  5.20:  Mesh  plot  after  deformation 


CHAPTER  6 

EXPERIMENTAL  PROGRAM 


The  first  section  of  this  chapter  will  discuss  the 
method  used  in  the  experimental  part  of  this  study  as  well 
as  the  formula  used  and  its  limitation.  The  other  section 
will  present  some  details  on  the  material,  type,  and  size  of 
specimen  used  in  addition  to  manufacturing  and  machining 
processes.  Finally,  a description  of  the  fixture  and 
mounting  procedure  will  be  given. 

6.1  Theoretical  Background 
The  experimental  program  was  designed  to  provide  some 
data  which  can  be  used  in  evaluating  the  strain  energy 
release  rate  G based  on  compliance.  This  method  is  also 
known  as  the  K-calibration  technique.  In  particular,  the 
program  was  intended  to  compare  an  experimental  evaluation 
of  G with  an  analytical  calculation  of  G based  on  the  stress 
intensity  factors  which  were  calculated  from  the  stress- 
hybrid  singularity  super  finite  elem.ent  (see  Section  4.6). 
This  comparison  leads  to  full  confidence  in  the  principles 
discussed  previously. 

This  compliance  calibration  technique  is  widely  used 
in  determining  the  fracture  toughness  of  isotropic  materials. 
This  method  is  not  valid  for  orthotropic  materials  when  the 
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crack  does  not  advance  along  the  original  crack  direction. 
This  is  because  the  compliance  calibration  technique  is 
based  on  the  principle  of  linear  fracture  mechanics,  which 
assumes  in  advance  that  the  crack  grows  along  its  original 
direction.  Because  of  the  above  restriction,  this  experi- 
mental verification  was  carried  out  for  a special  case  of 
two  S3rmmetric  cracks  emanating  from  a circular  hole  in  a 
direction  parallel  to  the  fiber  direction.  The  fiber  direc- 
tion was  chosen  to  be  perpendicular  to  the  loading  direc- 
tion in  order  to  produce  a single  mode  fracture  (opening 
mode) , so  that  one  could  measure  the  crack  opening  displace- 
ment (COD)  in  an  accurate  and  more  direct  way  by  using  the 
COD-extensometer . 

The  theoretical  expression  for  the  critical  strain 
energy  release  rate  using  the  compliance  calibration  method 
is  given  by  the  relation  [35] 


dX 


(6.1) 


where  (a-d/2)  is  the  original  crack  length  and  X is  the 
compliance  of  the  specimen,  which  is  defined  as  the  inverse 
of  the  slope  of  the  P-6  curve.  A polynomial  curve  was 
fitted  to  the  X vs.  a/w  data  and  the  slope  dX/da  was 
determined  as  a function  of  a/w.  The  energy  release  rate 
was  then  obtained  from  Equation  (6.1). 


110 


6.2  Experimental  Procedures 
6.2.1  Specimen  Preparation 

All  specimens  were  cut  from  two  rectangular  30.5  cm 
X 30.5  cm  plates  of  xmidirectional  composite  materials.  This 
composite  material  was  manufactured  from  glass-epoxy  which 
was  supplied  in  prepreg  Scotchply-1003  tape  rolls  by  the  3M 
Company  (Minnestoa  Mining  and  Manufacturing  Company).  The 
two  rectangular  plates  were  laid  up  in  one  direction  using 
15  plies  to  form  one  lamina  each.  They  were  pressurized 
and  cured  in  an  autoclave  using  the  prepreg  supplier's 
recommendation  for  the  curing  cycle  procedure.  For  more 
details  on  manufacturing  processes  and  a description  of 
experimental  equipment  see  Collazo's  study  in  reference 
[34].  The  elastic  constants  for  unidirectional  Scotch  ply 
glass-epoxy  are  the  following  [34]: 

= 40.0  GPa  (5.8  x 10®  psi) 

^2  “ 8.27  GPa  (1.2  x 10®  psi) 

^,2  = GPa  (0.6  X 10®  psi) 

V = 0.26 
1 2 

Six  identical  specimens  were  cut  and  machined  from  each 
plate  to  produce  a specimen  size  of  300  mm  x 42  mm.  Each 
specimen  was  weakened  by  drilling  a 10.66  mm  diameter  hole 
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in  the  center  of  the  specimen.  Two  S3nnmetric  cracks  were 
initiated  at  the  center  of  the  specimen  on  the  hole 
bovindary  using  a jeweler's  saw.  Some  specimens  suffered 
some  degree  of  delamination  near  the  hole  on  the  crack  side 
of  the  specimen  where  the  drill  exited,  while  the  center 
cracks  all  exhibited  some  degree  of  deviation  from  a 
straight  line  due  to  wandering  of  the  saw  blade. 

All  twelve  specimens  cut  from  the  two  plates  were 
identical  in  all  dimensions  and  fiber  orientation  except  for 
the  ratio  of  the  crack  length  to  specimen  width,  which  was 
nominally  0.3,  0.4,  0.5,  0.6,  0.7,  and  0.8  for  the  six 
specimens,  respectively.  For  illustration  of  all  specimen 
dimensions  see  Figure  6.1  The  other  six  specimens  were 
prepared  with  the  same  crack  length  to  specimen  width 
ratios;  thus,  there  were  two  identical  specimens  from  each 
type.  All  specimens  had  an  average  thickness  of  about 
3.25  mm  and  the  fiber  volume  fractions  for  the  plates  were 
determined  to  be  approximately  50  percent. 

6.2.2  Test  and  Procedures 

All  specimens  were  ramp  loaded  by  friction  grips 
(see  Figure  6.2)  and  tested  at  room  temperature  and  a 
constant  extension  rate  of  0.02  mm  per  second  to  the  maxi- 
mum load  Each  test  was  stopped  shortly  after 

reaching  the  maximum  load,  which  corresponded  to  the 
beginning  of  crack  propagation.  The  COD  versus  load  was 
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Fiber  direction 


w 21.0  mm;  L - 150  mm;  d = 10.66  mm; 

a/w  = 0.3,  0.4,  0.5,  0.6,  0.7,  0.8 


Figure  6.1;  A test  specimen  with  all  dimensions 
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Figure  6.2:  A specimen  mounted  by  friction  grips 
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recorded  for  all  specimens  by  using  an  extensometer  (MTS 
Model  632.02C-20),  which  was  moimted  on  the  center  of  the 
specimen  by  a special  fixture  (see  Figure  6.3).  This 
fixture  was  designed  to  attach  the  extensometer  on  the 
specimen  and  to  allow  it  to  be  aligned  centrally  on  the 
vertical  axes  of  the  specimen  for  better  recording  of  COD. 
Figure  6.4  shows  the  tension  specimen  with  extensometer 
attached,  momted  on  the  Material  Testing  System  (MTS) 
testing  machine.  The  applied  load  and  the  corresponding 
displacement  of  the  extensometer  were  recorded  continu- 
ously on  an  x.y-plotter  (Omnigraphic  2000;  Houston  Instru- 
ments) . A typical  load  versus  COD  curve  for  the  specimen 
having  a/w  — 0.5  is  shown  in  Figure  6.5.  The  load  which 

causes  incipient  crack  growth  is  defined  as  P and  the 

max 

corresponding  COD  is  the  maximum  crack  opening  displacement. 

The  inverse  of  the  slope  of  the  linear  portion  of  the  load 

versus  COD  is  known  as  the  compliance,  X,  of  the  load 

versus  COD  curve.  For  each  specimen,  data  of  P and  A 

max 

were  recorded. 
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Hex  nut : 8 required 


Flat  head 
screw: 

8 required 


Specimen 


(b)  • 


(c) 


(a)  MTS  extensometer  Model  632.02C-20 

(b)  Fixture  assembly 

(c)  Front  view  of  the  fixture  inside  the  hole 


Figure  6.3:  Extensometer  and  fixture  assembly 
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Figure  6.4:  Photographic  picture  of  MTS  testing  machine 
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Figure  6.5:  A typical  load  vs.  crack  opening  displacement 


CHAPTER  7 

EXPERIMENTAL  RESULTS,  CONCLUSIONS, 

AND  RECOr-IMENDATIONS 

This  chapter  will  present  the  test  results  of  glass/ 
epoxy  specimens  to  determine  the  critical  strain  energy 
release  rate.  These  results  will  then  be  compared  to  the 
numerical  results  obtained  from  the  finite  element  computer 
program.  Stress  intensity  factors  and  the  energy  release 
rate  will  be  obtained  for  the  case  of  a plate  with  two 
symmetrical  cracks  emanating  from  a circular  hole  in  uni- 
directional reinforced-f iber  composite  materials. 

Finally,  an  overall  conclusion  for  the  two  main 
parts  of  this  work  (stresses  around  an  opening  with  and 
without  cracks)  will  be  given.  The  last  section  will  be 
mainly  for  stating  some  suggestions  and  recommendations 
for  future  work  related  to  this  topic. 

7.1  Experimental  Results 

All  twelve  test  specimens  were  loaded  on  the  MTS 
testing  machine  to  failure.  Figure  7.1  shows  a typical 
photograph  of  a specimen  and  all  other  attachments.  A 
typical  record  of  all  specimens,  including  their  initial 
crack-length  to  width  ratio  (a/w) , the  compliance  X, 
the  derivative  dX/d(a/w),  and  the  critical  energy  release 
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Figure  7.1:  A typical  photograph  illustrating  a 

specimen,  fixture,  and  extensometer  all 
mounted  on  MTS  testing  machine 
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rate  per  unit  area  of  crack  surface  growth  (N-M/m^ 

= N/m),  all  of  which  have  been  tabulated  in  Table  7.1, 

A set  of  cubic  polynomials  (quasi-Hermite  spline) 
in  a form  of  mathematical  library  subroutine  [36]  was  used 
to  interpolate  a set  of  experimental  points  from  the  com- 
pliance vs.  initial  a/w  data.  This  quasi-Hermite  spline 
produces  a continuous  first  derivative  of  the  compliance 
dA/d(a/w).  This  first  derivative  was  evaluated  at  each 
initial  crack/width  ratio.  Then,  using  Equation  (6.1), 
the  critical  energy  release  rate  at  each  initial  crack/width 
ratio,  a/w,  can  be  obtained  as  in  Table  7.1. 

Data  obtained  from  the  experimental  evaluation  of 
compliance  were  scattered  above  and  below  the  curve  as  is 
indicated  in  Figure  7.2.  Therefore,  the  average  value  of 
these  data  was  used  to  fit  the  compliance  curve,  even 
though  this  average  seems  not  to  be  accurate  enough,  be- 
cause it  was  obtained  from  only  two  specimens.  Generally, 
these  experimental  results  give  us  some  confidence  about 
the  predicted  value  of  G^^,  from  the  analytical  method. 

The  experimental  value  for  the  critical  energy  release  rate 
for  glass/epoxy  was  determined,  for  the  case  of  a/w  = 

0,5,  to  be  18.3  N/m  where  its  calculated  value  is  39.90 
N/m.  This  difference  between  analytical  and  experimental 
evaluation  was  expected  for  three  reasons.  First  is  the 
small  number  of  specimens  used.  Second  is  the  difficulty 


Table  7.1:  Experimental  data  and  results  for  obtaining  energy  release  rates 
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Crack  length  ratio  (a/w) 


Figure  7.2:  Compliance  A vs.  a/w  plot 
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of  producing  identical  specimens  due  to  manufacturing  and 
machining  all  specimens  manually  in  the  laboratory.  The 
third  reason  is  the  difficulty  of  producing  a sharp  crack 
tip  as  in  the  theoretical  analysis. 

7.2  Analytical  Results 

To  produce  some  numerical  results  we  chose  the  case 
of  two  initial  symmetric  cracks  perpendicular  to  the  load- 
ing direction  and  emanating  from  a circular  hole  with 
ratio  a/w  = 0.5,  for  five  different  fiber  orientations 
with  respect  to  the  loading  axis.  In  the  case  where  the 
fiber  direction  angle  is  0°  or  90° , the  calculation  of  the 
stress  intensity  factors  and  the  energy  release  rate  were 
carried  out  for  one  quarter  of  the  specimens.  This  quarter 
was  divided  into  29  quadrilateral  isoparametric  elements, 
one  being  the  singular-super  element,  where  the  stress 
singularity  near  the  crack  tip  has  been  treated  by  using 
the  hybrid  stress  model  discussed  in  Chapter  4.  The  total 
number  of  degrees  of  freedom  was  81;  Figure  7.3  shows  a 
computer  plot  mesh  for  this  case.  Calculations  have  been 
carried  out  for  glass/epoxy  and  graphite/epoxy  composite 
materials . 

For  a case  where  the  fiber  direction  makes  some 
angle  with  the  loading  direction,  we  considered  the  full 
specimen  mesh,  because  there  was  no  material  or  geometric 
symmetry.  The  specimen  was  divided  into  114  quadrilateral 
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isoparametric  elements,  two  of  which  were  singular-super 
elements,  one  for  each  crack  tip.  The  total  number  of 
degrees  of  freedom  was  322.  Computer  mesh  plots  for  the 
full  specimen  before  and  after  deformation  are  shown  in 
Figures  7.4  and  7.5,  respectively. 

All  results  obtained  by  the  finite  element  computer 
program  are  tabulated  in  Tables  7.1  and  7.2  for  glass/ 
epoxy  and  graphite/epoxy  composite  materials,  respectively. 
Figure  7.6  shows  the  result  of  Kj.  and  with  different 
fiber  orientations  for  glass/epoxy  composite  material,  as 
given  by  Equations  (4.41)  and  (4.42)  of  Section  4.5. 

Figure  7.7  shows  the  results  of  calculation  of  strain 
energy  release  rate  for  glass/epoxy  having  different  fiber 
orientations,  as  given  by  Equation  (4.44)  of  Section  4.6. 

7.3  Overall  Conclusions  and  Remarks 
7.3.1  Stress  Concentration  Factors 

Based  on  the  analytical  solution  discussed  in 
Chapter  2,  the  results  shown  in  Figures  5.1  through  5.18 
and  Tables  5.1  through  5.4,  we  make  the  following  remarks 
regarding  the  stress  concentration  near  holes  in  unidirec- 
tional fiber-reinforced  composite  materials : 

!•  The  maximum  value  of  stress  concentration 
depends  on  the  stiffness  of  the  composite. 
The  location  of  this  maximum  value  depends 
on  the  fiber  orientation,  fiber  material, 
and  loading  direction. 


2. 


126 


i i 

1 

1 

1 

1 

li 1 

II 

Figure  7.4:  Full  specimen  mesh  before  deformation 

(nonsymmetric  case),  a/w  = 0.5 


L27 


Figure  7.5:  Full  specimen  mesh  after  deformation 

(nonsymmetric  case),  a/w  = 0.5 
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Table  7.2:  Analytical  results  for  K^. , K^j.,  and  G 


Fiber  Orientation 
0 (degree) 

Kj  xlO-® 
(MN/m^)  i/m 

Kjj  xlO-® 

(MN/m^) 

G 

(N/m) 

0 

9.07 

0.00 

261.20 

20 

3.62 

7.30 

156.00 

45 

1.03 

6.35 

84.27 

60 

3.24 

4.04 

44.94 

90 

5.26 

0.00 

39.90 

Note:  These  results  obtained  from  the  finite  element  com- 

puter program  for  the  case  of  unidirectional  glass/ 
epoxy  composite  material 

Table  7.3:  Analytical  results 

for  Kj,  Kj.j,  and  G 

Fiber  Orientation 
0 (degree) 

Kj  X 10'* 

(MN/m^)v^ 

Kjj  xl0‘^ 

(MN/m^)\/m 

G 

(N/m) 

0 

9.57 

0.00 

194.90 

90 

3.92 

0.00 

9.22 

Note:  These  results  obtained  from  the  finite  element  com- 

puter program  for  the  case  of  unidirectional  graphite/ 
epoxy  (symmetric  case) 


Stress  intensity  factors  (K 


129 


Figure  7.6:  Stress  intensity  factors  vs.  fiber  orientation 

for  two  symmetric  cracks  emanating  from  a circular 
hole  in  unidirectional  glass/epoxy  (a/w  = 0.5) 


Energy  release  rate  (G)  N/m 
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Fiber  orientation  (©) 


Figure  7.7:  Energy  release  rate  vs.  fiber  orientation 

for  glass/epoxy  (a/w  = 0.5) 
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3.  The  maximvim  stress  can  be  reduced  by 

changing  the  fiber  direction  with  respect 
to  loading  direction. 

A.  The  shape  and  size  of  the  hole  also  has  a 
large  effect  on  stresses.  For  example,  a 
hole  of  a circular  shape  causes  less 
stress  than  an  elliptical  shape  with  the 
major  axis  perpendicular  to  the  loading 
axis . 

5.  When  the  value  of  the  aspect  ratio  a/b 
for  the  elliptical  hole  increases,  the 
location  of  maximum  stress-concentration 
becomes  more  independent  from  both  fiber 
orientation  and  fiber  material  for  the  same 
loading  direction. 

6.  As  the  aspect  ratio  (a/b)  increases,  the 
maximiim  stress-concentration  factor 
increases  linearly  for  glass/epoxy  and 
graphite/epoxy  composite  materials.  The 
same  is  also  true  for  the  isotropic  case. 

7.3.2  Stress  Intensity  Factors 

The  stress  intensity  factor  for  the  opening  mode 
with  the  fiber  direction  parallel  to  the  loading  direction 
has  a higher  value  than  the  value  for  the  case  with  the 
fiber  direction  perpendicular  to  the  loading  direction.  On 
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the  other  hand,  uncracked  material  in  the  first  case  with 
the  fibers  parallel  to  the  loading  direction,  is  considered 
to  be  much  stronger  than  in  the  other  case. 

The  same  argument  holds  for  the  energy  release 
rate,  where  we  have  a higher  value  for  the  case  of  fiber 
direction  being  parallel  to  the  loading  direction.  In 
comparing  the  single  mode  fracture  behavior  between  the 
two  different  composite  materials  studied,  we  find  that  the 
graphite/epoxy  composite  gives  a lower  stress  intensity 
factor  and  a lower  energy  release  rate  than  glass/epoxy. 
This  means  that  a composite  material  made  of  graphite/ 
epoxy  offers  more  fracture  resistance  than  a composite 
material  made  of  glass/epoxy. 

In  the  mixed-mode  fracture  behavior,  the  fracture 
failure  takes  place  because  of  a combination  of  two  stress 
intensity  factors,  and  In  the  analytical  results 

of  Figure* 7.6,  we  find  that  the  stress  intensity  factor 
of  the  shear  mode  decreases  as  the  fiber  direction  angle  0 
increases.  But  increasing  this  angle  reduces  the  strength 
of  the  material  in  the  loading  direction.  So,  it  seems 
impossible  to  optimize  and  at  the  same  time.  Obvi- 
ously, there  is  a trade-off  between  K^.  and  depending 

upon  the  loading  direction  and  design  requirements. 
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7.4  Recommendations  for  Future  Study 

There  are  two  distinct  areas  in  which  future  research 
can  be  carried  out.  The  first  deals  with  the  failure  cri- 
teria and  the  location  of  failure  stresses  due  to  the  high 
stress-concentrations  near  holes  in  composite  materials. 

In  addition,  a study  of  the  effects  of  the  stacking  sequence, 
and  the  type  of  material  and  loading  condition  on  the 
fracture  properties  of  a laminate,  will  make  a valuable 
contribution  to  engineering  design. 

The  second  general  area  of  future  research  is  to 
study  the  stress  intensity  factors  using  different  finite 
element  models  based  on  one  of  the  mixed  energy  approaches 
discussed  in  Chapter  3 and  then  compare  its  results  with 
the  results  obtained  from  the  approach  studied  here. 

Another  possible  continuation  of  this  work  is  to  extend 
the  technique  adopted  in  this  dissertation  to  three- 
dimensional  elasticity  and  then  solve  the  fracture  mechanics 
problem.  A special  treatment  has  to  be  added  into  the 
code  in  order  to  seek  a solution  when  we  have  pairwise 
equal  roots,  such  as  were  obtained  by  solving  the  charac- 
teristic equation  (2.11)  of  Chapter  2. 

One  final  recommendation  is  to  study  the  order  of 
the  stress  singularity  near  the  crack  tip  for  different 
composites  and  its  effect  on  the  evaluation  of  stress 
intensity  factors  near  the  crack  tip. 
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APPENDIX  A 

ELASTIC  COEFFICIENTS  FOR 
GEIJERAL  ORTHOTROPIC  MATERIALS 


The  coefficients  of  the  stress-strain  matrix  a. . of 

ij 

Equation  (2,4)  can  be  written  as  follows  [18,  37]: 


a = cos‘*0/E  + (1/G  - 2v  /E  ) sin^e  cos^e  + sin‘*e/E 

11  1 12  12  1 2 


a = - V (sin'^e  + cos‘*6)/E  + (1/E  + 1/E 

12  12  112 


- 1/G  ) sin^e  cos ^6 

1 2 


a = Sin‘*e/E  + (1/G  - 2v  /E  ) sin^e  cos^e  + cos‘*e/E 

22  1 12  12  1 2 


a = (2/E  + 2v  /E  - 1/G  ) Sine  cos^e 

1 e 1 121  12 


- (2/E  + 2v  /E  - 1/G  ) Sin^e  cos6 

2 121  12 


a = (2/E  + 2v  /E  - 1/G  ) sin^B  cosB 

26  1 121  12 


- (2/E  + 2v  /E  - 1/G  ) sine  Cos^B 

2 121  12 
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a = (4/E  + 4/E  + 8v  /E  - 2/G  ) sin^e  cos^e 

66  1 2 12  1 12 

+ (sin ‘*6  + cos‘*6)/G 

12 

where  9 is  the  angle  between  the  fiber  direction  and  the  x- 
axis  (see  Figure  2,1).  All  the  above  coefficients  are 
expressed  in  terms  of  just  four  independent  material  con- 
stants E , E , V , and  G 


APPENDIX  B 

THE  MATHEMATICAL  FOP>MULATION  OF  THE 
CONVENTIONAL  FINITE  ELE>iENT 
(DISPLACEMENT  MODEL) 

For  the  finite  element  displacement  model,  the 
governing  equilibrium  equations  can  be  obtained  by  minim*iz- 
ing  the  total  potential  energy  functional  of  the  system.  The 
total  potential  energy  functional  (ir)  can  be  expressed  as 
previously  stated  in  Equation  (3.2) 

TT  = I i E . . , . R . . e,  . dV 

p J 2 ijk£  ij  k£, 

V 

or  in  matrix  notation 

TTp  = f ^ {e}^[E]{£}  dV  - f {u}^{T}  dS  (B.2) 

V S 

a 

In  the  finite  element  displacement  method,  the  displacement 
is  assumed  to  have  unknown  values  only  at  the  nodal  points, 
so  that  the  variations  within  any  element  are  described  in 
terms  of  the  nodal  values  by  means  of  interpolation  function. 
Thus , 

{u}  = [L] {q®}  (B.3) 


T.  u.  dS  (B.l) 
s 

a 
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where  IL]  is  the  set  of  interpolation  functions  termed 


ments  of  the  element.  The  strains  within  the  elem.ent  can 
be  expressed  in  terms  of  element  nodal  displacement  as 


where  [B]  is  the  strain  matrix  generally  composed  of 
derivatives  of  the  shape  functions.  If  we  assume  that  the 
total  volume  V is  divided  into  M discrete  elements  and  the 
element  shape  functions  have  been  chosen  so  that  no  singu- 
larities exist  in  the  integrands  of  the  functional,  the 
total  functional  of  potential  energy  of  the  continuum  can 
be  expressed  as  the  sum  of  the  functional  contributions  of 
the  individual  element . Thus , 


where  tt^(u,v)  represents  the  potential  energy  functional  of 
an  element  e which  can  be  written 


shape  functions  and  {q®}  is  the  vector  of  nodal  displace- 


{£} = [BHq®} 


(B.4) 


M 

TI  (u,v)  = I 7r^(u,v) 
® e = l 


(B.5) 


(B.6) 
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We  know  that  at  equilibritun  the  potential  energy  of  the 
system  assumes  a minimiam  value.  Because  of  the  summation 
principle  expressed  in  Equation  (B.5),  we  may  carry  out  the 
minimization  process  element  by  element.  We  note  that 
the  potential  energy  of  the  discretized  system  assumes  its 
minimum  value  when  the  first  variation  of  the  functional 
vanishes,  that  is, 


M 

6tt  (u,v)  = I 671  (u,v)  = 0 (B.7) 

P e=l 


where 


6tt  (u,v)  = I 


. , 9u. 
1 = 1 1 


377  r 377 

' ''“i  87T  ° 

i = l 1 


(B.8) 


and  r is  the  niomber  of  element  nodes.  But  the  3u.  and  the 

1 

3v^  are  independent  variables;  hence,  we  must  have 


3T7  3T7 

e _ e 

3u.  3v. 

1 1 


i=  1 2 r 


(B.9) 


for  every  element  e of  the  system.  Thus,  after  performing 
this  minimization  process  on  Equatipn  (B.6),  it  becomes 


[B]  [E] [B] {q®}  dV  - 


[L]^{T}  ds  = 0 


or 


[K®]{q®}  - {F®}  = 0 


(B.IO) 


140 


where 


IK®] 


[B]^[E][B]  dv 


(B.ll) 


V 


e 


is  the  "element  stiffness  matrix,"  and 


{F®} 


[L]'^{T}  ds 


(B.12) 


S 


e 


are  the  equivalent  nodal  forces  for  the  element.  The  sum.- 
mation  of  the  terms  in  Equation  (B.IO)  over  all  the  elements 
results  in  a system  of  equilibrium  equations  for  the  entire 
continuum.  These  equations  are  then  solved  by  any  standard 
technique  to  yield  the  nodal  displacements  (q). 


APPENDIX  C 

STRESS  AND  DISPLACEMENT  EQUATIONS 
FOR  NEAR  CRACK  TIP  SOLUTION 

The  stress  field  near  the  crack  tip  in  an  anisotropic 
body,  as  was  stated  by  Sih,  Paris,  and  Irwin  [32] , can  be 
divided  into  modes  of  deformation,  but  the  degree  of  simpli- 
fication thus  achieved  is  less  than  for  the  isotropic  case. 
This  is  because  the  crack  surface  displacements  in  general 
anisotropic  media  depend  upon  the  directional  properties  of 
the  material  and  do  not  necessarily  occur  in  a planar  fashion. 
However,  the  general  state  of  stress  and  displacement  near  a 
crack  tip  can  still  be  considered  as  the  sum  of  three  indi- 
vidual boundary  problems,  namely,  those  corresponding  to 
inp lane -sjnnme trie  loading,  inplane- skews)mime trie  loading, 
and  antiplane  shear. 

The  stresses  and  displacments  equations  in  a small 
region  surrounding  the  crack  tip  in  plane  deformation  are 
obtained  from  the  same  reference  [32]  as  follows: 
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Symmetric  loading  (crack  opening,  mode  I): 
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c 

to 


Skew- symmetric  loading  (crack  sliding,  mode  II): 


143 


c 

CO 


CNJ 
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The  stress  distributions  and  the  displacement  field 
near  the  crack  tip  in  an  isotropic  body  can  be  expressed  in 
the  following  due  to  the  work  of  Irwin  [38]  and  Westergaard 
[39]  : 

1.  Symmetric  loading  (crack  opening,  mode  I): 


a = — cos  6/2(1  - sin  6/2  sin  36/2) 

^ /57 


o = — cos  6/2(1  + sin  6/2  sin  36/2)  (C.5) 

^ /2? 


T = — - COS  6/2  sin  6/2  cos  36/2 


and 


u = (1+v)[(2k  -1)  cos  6/2  - cos  36/2] 


2E 


(C.6) 


= ~/j  (1+v)[(2k+1)  sin  6/2  - sin  36/2] 


V 
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2.  Skew- symmetric  loading  (crack  sliding,  mode  II): 


o = - — sin  0/2  (2  + cos  6/2  cos  36/2) 
""  /2? 


a 

y 


— ^ sin  6/2  cos  6/2  cos  36/2 


T = -5-^-  cos  6/2  (1  - sin  6/2  sin  36/2) 


and 


u = (1  + v)[(2k+3)  sin  6/2+sin  "iQll] 

2E  V ^ 


V = - /y  (1  +v)  [(2k  - 3)  cos  6/2  - cos  36/2] 

2E  * ^ 


where  k = (3  - v)/(l  + v)  for  plane  stress. 
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DERIVATION  OF  EQUATION  (4.1) 


Starting  from  the  modified  complementary  energy 
functional  of  Chapter  3 and  then  adding  and  subtracting 
the  term 


to  Equation  (3.8)  with  multipliers  replaced  by  their 
physical  meaning  u^,  interelement  boundary  displacements, 
Equation  (3.8)  we  derive 


n 


n 


0, 


u 


n 


n 


+ 


(D.l) 


a, 


n 


let 
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9A 


Ti  dS  = 


dS  + 


Ti  Ui  dS  + 


n 


n 


u 


n 


Ti  Ui  dS 


(D.2) 


where  9A  = S + S + S is  the  entire  boundary  of  an  element 

n n 
A 

n 


u.  . u. 


Ti  - T, 


on  S 


u 


on  S 


n 


n 


(D.3) 


Also  let 


M 


= I 


(D.4) 


n=l 


where  tt  now  denotes  the  complementary  energy  functional  for 
one  element  A with  boundary  3A. 

Substitution  of  Equations  (D.2),  (D.3),  and  (D.4) 
into  Equation  (D.l)  yields 
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1 p 

2 ^ijk£  ^ij 


dA  - 


~ r _ ~ 

T.  u.  dS  + T.  u.  dS 

J 9A  ^ ^ •'S  ^ 

o 

(D.5) 


The  Euler  equations  of  the  functional  of 

Equation  (3.S)  of  Chapter  3 are 


1A8 


Cijkt  “M  ' I ("i,j  “j,l' 


(D.6a) 


o • • . = 0 

ij  ,3 


in  A 


(D.6b) 


T . = o . . V • 
1 ij  J 


on  9A 


(D.6c) 


Now,  with  the  help  of  Equation  (D.6)  and  the  divergence 
theorem,  let  us  take  a look  at  the  first  term  on  the  right 
hand  side  of  Equation  (D.5)  as  the  following 


1 


C . „ o . • Oi.  „ dA  = 


I ^ijk£  ‘^ij  ^k£ 


1 

I 


^ki  °ki 


“ij  “j.i>‘^ 


(D.7) 


But 


dA  = 


[(aij 


“i>.j 


- a 


ij  ,J 


u^]dA 


Now,  applying  the  divergence  theorem  to  the  first 
term  on  the  right  hand  side  of  the  above  equation,  and 
cancelling  the  second  term  due  to  Equation  (D.6b)  leads  to 

f 

a.,  u.  . dA  = o..  u.  V.  dS  (D.8) 

Ja  JgA  ^ J 


Inserting  Equation  (D.6c)  into  Equation  (D.8),  we  get 
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a. . u.  . dA  = T.  u.  dS  (D.9a) 

Ja  J 3A 

and  similarly 

o. . u.  . dA  = 1 T.  u.  dS  (D.9b) 

Ja  J.i  J 9A  ^ 

Insertion  of  Equations  (D.9a)  and  (D.9b)  into  Equation  (D.7) 
gives 


1 

J 


•A 


u- 


dS 


(D.IO) 


By  replacing  the  first  term  of  Equation  (D.5)  with  its 
expression  from  Equation  (D.IO)  we  obtain  the  following: 


dA 


T.  u.  dS 


dS  + 


dS 


but  the  term 


T.  u.  dS  = 0 
Js  11 
o 


because  of  the  free  traction  on  the  crack  surfaces. 


APPENDIX  E 

FINITE  ELEMENT  COMPUTER  PROGRAM 


On  the  following  pages  appears  a complete  listing  of 
the  finite  elements  computer  program  which  has  been  developed 
and  used  in  this  research.  This  program  was  written  in  a 
standard  FORTRAN  language  and  is  applicable  to  FORTRAN 
compiler  G and  all  other  modified  versions  of  the  FORTRAN 
compilers . 

The  listing  of  the  matrix  inversion  and  the  plot 
subroutines  are  not  included  here.  The  matrix  inversion 
subroutine  was  taken  from  the  Standard  International  Mathe- 
matical and  Statistical  Libraries  (IMSL) , while  the  plot 
subroutines  were  taken  from  the  Gould  plot.  Both  of  these 
subroutines  were  installed  by  the  Northeast  Regional  Data 
Center  (NRDC) , State  University  System  of  Florida. 

This  finite  elements  computer  code  handles  the  two- 
dimensional  elasticity  problem  (plane  stress)  in  fracture 
mechanics.  This  code  calculates  stress  intensity  factors 
(Kj.,Kj^j)  near  the  crack  tip  in  the  presence  of  a single  or 
mixed  mode  fracture  for  either  isotropic  or  anisotropic 
material . In  addition  it  calculates  the  energy  release  rate 
for  the  case  of  rectilinear  anisotropic  materials. 


150 


151 


IHPLICIT  REAL«8(A-HiO-Z) 

COMPLEX  ZHU1>ZMU2>RMU1>RMU2 

REAU4  Y0UNGlfY0UNG2fP0IS12.G12»THICK»SP(6)fXTIP(2)»YTIP(2) 

% »R0TR(4)>R0TI(4) 

DIMENSION  COORD(72i2)»LNODS(46i9)fASDIS(144)»NOFIX(  2)i 
iIFPRE(  2»2)fPRESC(  2»2)  ,EL0A0(46»18) »NDFR0(72) »BCRT(2f2il8) 
t »CRD<72»2)fSIF(4) 

COMMON/WORK/SHAP ( 4 ) i DERI V ( 2 » 4 ) » BMATX  < 3 » 8 ) » POSGP ( 3 ) 
t iUEIGP(3)»QBHAT(3f8)»GPC0D(2>9)fSMATX(3f8>9) 

COMMON/MATR I L/YOUNGl » Y0UNG2 » POI S 1 2 » G 1 2 1 1 ANGL » I STRPC » TH I CK 
COMMON/SUPER/ ISUPR 1 » I SUPR2  f KE Y 1 » KE Y2 1 XT I P » YT I P 
C 

REWIND  1 
REWIND  2 

PI=4,0»DATAN(1.0D0) 

READ<5f900)NP0IN»NELEHfNVFIXf  NGAUSiNN0DE»ISUPRl»ISUPR2fKEYl»KEY2 
WRITE(6»905)NP0IN»NELEMiNVFIXiNGAUSiNN0DE»ISUPRl»ISUPR2fKEYl»KEY2 
NT0TV=2*NP0IN 

Ct«»<«Re3d  all  the  eleaent  nodal  coordinate  (Xfv)  

DO  5 I<1>NP0IN>6 

12=1+5 

READ(5.930)(<C00RD(KK»JJ)»JJ=l»2)fKK=I»I2) 

5 CONTINUE 

DO  3 I=lfNELEM 
DO  3 J=liNNODE 
3 LN0DS(I»J)=0 

DO  10  I=lfNELEM»4 
K=I+3 

IF(I.EQ»45)  K=46 

READ(5»925)((LN0DS(KK»J)f J=1.4).KK=I»K) 

10  CONTINUE 

RE AD ( 5 » 926 ) ( LNODS ( I SUPR 1 » I ) » I =5 1 9 ) 
READ(5»926HLN0DS(ISUPR2»I)tI=5»9) 

WRITE(6»915) 

DO  20  IELEM=1»NELEM 

20  WRITE(6»920)  IELEM»  (LNODS(IELEM»INODE)f IN0DE=1.NN0DE) 
WRITE(6>940) 
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URITE(6>945) 

C Read  and  write  the  fixed  values  (Restranined  nodes) 

DO  60  IVFIX=lfNVFIX 

READ(5»950)  N0FIX<IVFIX)»(IFPRE(IVFIX»ID0FN)»ID0FN=1»2) 
f » ( PRESC ( I VFIX » IDOFN ) , ID0FN=1 1 2 ) 

60  «RITE(6f951)  NOFIX(IVFIX) » (IFPREdVFIXf  IDOFN)  * ID0FN=1 .2) 
f »(PRESC(IVFIX»ID0FN)»ID0FN=1»2) 

C<««««Set  UP  the  saeple  GAUSS  point  positions  and  UeiShtinS  Factors 
C for  nuarical  integration. 

POSGP  < 1 ) =-0 . 77459666924 1 4B3D0 

P0SCP(2)=0,0D0 

P0SGP(3)=-P0SGP(1) 

UEIDP( 1 )=0 .555555555555556D0 
VE IGP  ( 2 ) :=0 . 888888888888889D0 
HEIGP<3)=WEIGP(1) 

C 

READ(5»960)  YOUNGlfYOUNG2»POIS12fG12»IANGL»ISTRPC»THICK 
WRITE(6i965)Y0UNGl»Y0UNG2»P0IS12»G12»IANGLfISTRPC»THICK 
C 

C Read  Crack  Tip  Diiiensions 

REAIU5f800)  XTIP(l)fYTIP(l)»XTIP(2)»YTIP(2) 
WRITE(6f805)XTIP(l)»YTIP(l)»XTIP(2)»YTIP(2) 

800  FORMAT(4E10.0) 

805  F0RHAT(3Xf'XTIP(l)='»F7.3»4X.'YTIP(l)='fF7.3.4X» 

J 'XTIP(2)='fF7.3f4X»'YTIP(2)='»F7.3»//> 

C 

CALL  CHECK ( NELEM » NVF I X » NOF I X » LNODS » I FPRE » COORD i NPO IN » 

» NFrSNfNDFRO) 

C NEXT  CREAT  THE  ELEMENT  STIFFNESS  FILE 

CALL  STIFPS(NELEM*NPOIN»NGAUS*COORDfLNODSfBCRTfSP»ROTR»ROTI) 

C COMPUTE  LOADS 

CALL  LOADPS ( NELEM . NPOIN » NGAUS  f COORD  > LNODS  f ELOAD ) 

C SOLVE  THE  RESULTING  EQUATIONS  BY  THE  FRONTAL  SOLVER.... 

C 

CALL  FRONT (NELEM»  NPOINf  NVFIX  > NTOTV » NGAUS  t COORD » LNODS*  ELOAD» 

* ASD IS  * NOFIX  * I FPRE  * PRESC  * BCRT  * MFRON  * CRD  * S I F ) 
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C CALCULATE  ENERGY  RELEASE  RATE  (ANISOTROPIC  CASE) 

IF(ISTRPC.EQ.l)  GO  TO  37 
ZHU 1 =CHPLX ( ROTR ( 2 ) » ROT I ( 2 ) ) 

ZHU2=CMPLX(R0TR(4) fR0TI(4) ) 

RMU1=1.0/ZMU1 

RMU2=1.0/ZMU2 

TEMJ1=SIF(1)*AIMAG(RHU2+RMU1)+SIF(2)»(REAL(RMUD* 

% AIHAG ( RKU2) (REAL (RNU2  >*AIHAG ( RHUl ) ) 

ENG1=-0.5»PI*SIF(1)»SP<3)»TEHJ1 

TENJ2=SIF(2)»AIKAG<ZHUHZMU2)+SIF(1)*(REAL(ZMU1)*AIHAG(ZHU2)+ 
i REAL(ZHU2)»AIHAG(ZNU1)) 

ENG2=0 . 5*PI »SIF (2 ) >SP ( 1 ) »TEM J2 
TOTAL=ENGHENG2 
URITE(6>970)  ENGl »ENG2> TOTAL 
37  CONTINUE 

970  F0RHAT(//»3Xi'Jl='»E12.4.6X.'J2='.E12.4»//. 

t 10X»'ENERGY  RELEASE  RATE=' »E12.4 ) 

900  F0RMAT(9I5) 

905  F0RMAT(//.8H  NPOIN  =»I4»4Xf8H  NELEH  =tI4f4X»8H  NVFIX  =iI4 

I f 4X f 8H  NGAUS  = » 14 1 4X » ' NNODE= ' » 1 4 » 4X » ' ISUPER 1 = '»I3»4X»'I SUPER2= ' t 
» I3»4Xf 'KEYl='fI3»4Xf 'KEY2='»I3f//> 

915  FORMAT <//»5X»8H  ELEMENT »6X»13H  NODE  NUMBERS) 

920  F0RMAT(lX»I5i6X»9I5) 

925  F0RMAT(16I4) 

926  FORMAT (514) 

930  F0RMAT(12D6.0) 

940  F0RMAT(//»5X»17H  RESTRAINED  NODES) 

945  F0RMAT(5H  N0DE»1X.4HC0DE»6X»12HFIXED  VALUES) 

950  F0RMAT(1X»I4»3X»2I1»2D8.5) 

951  F0RMAT(1X.I4.3X.2I1»2D14.5) 

960  FORMAT(4E7,0»2I7»E7.0) 

965  F0RMAT(/»3X» 'Y0UNG1=' »E10,3»5Xf 'Y0UNG2=' .E10.3»5Xr 'P0IS12=' »F7.5» 
i6Xf'G12='»E10.4f6X»'ANGLE='»I4»4Xf 'ISTRPC='.Il»4X»'THICK='fF5.2f/) 
STOP 
END 
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C 


SUBROUT I NE  CHECK ( NELEH » N VF I X » NOF IX » LNODS » I FPRE  t COORD i NFO I N » 
MFRON.NDFRO) 


IMPLICIT  REAL»8<A-HfO-Z) 

DIMENSION  NDFRO(NELEM)iNOFIX(NUFIX)fLNODS(NELEH>9)f 
1 IFPRE(NVFIXf2)»NER0R(8)fC00RD(NP0INi2) 

REALtA  XTIP(2)fYTIP(2) 

C0MM0N/SUPER/ISUPR1»ISUPR2.KEY1.KEY2»XTIP,YTIP 
DO  10  I^lfNELEH 
10  NDFR0(I)=0 
DO  15  I=li8 
15  NEROR(I)sO 

C CHECK  FOR  ANY  REPETITION  OF  A NODE  NUMBER  WITHIN  AN  ELEMENT 

DO  140  IP0IN=1»NP0IN 
KSTAR=0 

DO  90  I=1»NELEM 

KZER0=0 

NN0DE=4 

IF ( I . EQ . ISUPRl . AND . KEYl . EQ . 1 )NN0DE=5 
IF ( I . EQ . ISUPRl . AND . KEYl . EQ , 2 ) NN0DE=9 
IF(I.EQ.ISUPR2)NN0DE=9 
DO  90  IN0DE*1*NN0DE 
IF(LNODS(I»INODE).NE.IPOIN)  60  TO  90 
KZER0=KZER0+1 

IF ( KZERO . GT . 1 ) NEROR ( 1 ) =NEROR ( 1 ) +1 
IF<KSTAR.NE.O)  GO  TO  80 
KSTAR=I 

C CALCULATE  INCREASE  OR  DECREASE  IN  FRONTUIDTH  AT  EACH  ELEMENT  STAGE. 

NDFR0(I)=NDFR0(I)+2 
80  CONTINUE 


C AND  CHANGE  THE  SIGN  OF  THE  LAST  APPREANCE  OF  EACH  NODE 

KLAST=I 
NLAST=INODE 
90  CONTINUE 

IF(KLAST.EQ.O)  GO  TO  110 

IF ( KLAST .LT . NELEM ) NDFRO ( KLAST+1 ) «NDFRO ( KLAST4 1 ) -2 
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LMODS (KLAST  t NLAST ) «-IP0IN 
60  TO  HO 

C CHECK  THAT  COORDINATES  FOR  AN  UNUSED  NODE  HAVE  NOT  BEEN  SPECIFIED. 

110  yRITE(6>900)  IPOIN 

900  F0RHAT(/f3Xfl5H  CHECK  WHY  M0DEfI4»HH  NEVER  APPEARS) 
NER0R(2)sNER0R(2)H 
SIGKA=0.0 

DO  120  IDIHE=li2 

120  SIGHA=SIGHA4DABS(C00RD( IPOINt IDIHE) ) 
IF(SIGHA,NE.0.0)NER0R(3)=NER0R(3)+1 

C CHECK  THAT  AN  UNUSED  NODE  NUMBER  IS  NOT  A RESTRAINED  NODE 

DO  130  IVFIX=lfNVFIX 

130  IF(N0FIX<IVFIX).E0.IP0IN)NER0R(4)=NER0R(4)+1 
HO  CONTINUE 

C CALCULATE  THE  LARGEST  FRONTHIDTH. 

HFR0N=0 

MFR0N=0 

DO  150  HliNELEH 
NFRON=NFRON+NDFRO(I) 

150  IF(NFRON.GT.MFRON)  HFRON=NFRON 
WRITE (6f 905)  MFRON 

905  FORMAT (3X»/»29H  MAX  FRONTMIDTH  ENCOUNTERED  =»I5»/) 

C CONTINUE  CHECKING  THE  DATA  FOR  THE  FIXED  VALUES. 

DO  171  IVFIX=1.NVFIX 

IF(NOFIX(IVFIX).LE.O.OR.NOFIX(IVFIX).GT.NPOIN)NEROR(6) 
f =NER0R(6)F1 

K0UNT*0 

DO  160  ID0FN=1>2 

160  IF(lFPRE(IVFIXflD0FN).6T.0)  K0UNT=1 
IF ( KOUNT . EQ . 0 ) NEROR ( 7 ) =NEROR ( 7 ) +1 
KVFIX=IVFIX-1 
IF(KVFIX.E0.0)G0  TO  171 
DO  170  JVFIX=1»KVFIX 

170  IF(IVFIX.NE.1.AND.N0FIX(IVFIX).EQ.N0FIX(JVFIX))NER0R(8) 

f =NER0R(8)+1 

171  CONTINUE 
KER0R=0 
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DO  ISO  IER0R=lf8 
IF(NERORdEROR).EQ.O)  GO  TO  180 
KER0R=1 

WRITE(6f910)IER0R>NER0R(IER0R) 

910  FORMAT (//»3X»'m«  DIAGNOSIS  BY  CHECK»  ERROR  '»I3» 

» 18H  ASSOCIATED  NUMBER >15) 

180  CONTINUE 

IF(KEROR.NE.O)  GO  TO  200 

C..... RETURN  ALL  NODAL  CONNECTION  NUMBERS  TO  POSTIVE  VALUES 
DO  190  I=ltNELEM 
NN0DE=4 

IF  < I . EO . ISUPR 1 . AND . KE Y1 . EQ . 1 ) NN0DE=5 
IF ( I . EQ . I SUPR 1 . AND . KE Y 1 , EQ . 2 ) NN0DE=9 
IF(I.EQ.ISUPR2)NN0DE=9 
DO  190  IN0DE=1»NN0DE 
1 90  LNODS ( I » INODE ) =I ABS ( LNODS ( I » INODE ) ) 

RETURN 
200  STOP 
END 

SUBROUTINE  SHAPE (S.T) 

IMPLICIT  REALtS(A-HiO-Z) 

C 

c CALCULATES  SHAPE  FUNCTIONS  <SHAP)  AND  THEIR  DERIVATIVES  (DERIV) 

C FOR  2D  ISOPARAMATRIC  QUADRILATERAL  ELEMENT  (4  NODES). 

C 

DIMENSION  SI(4)>TI(4) 

DATA  SIfTI/-.5D0».5D0».5D0»-.5D0i-.5D0»-.5D0r.5D0».5D0/ 
COHMON/UORK/SHAP ( 4 ) > DERI V ( 2 > 4 ) » BMATX ( 3 » 8 ) » POSGP ( 3 ) 
f rWEIGP(3)>QBMAT<3>8)>GPC0D(2f9)>SHATX(3>8i9) 

C 

C COMPUTE  SHAPE  FUNCTIONS  AND  DERIVATIVE  IN  NATURAL  COORDINATES.... 

DO  10  1=1*4 

SHAP(I)=(0.5D0+SI(I)»S)t(0.5D0+TKI)»T) 

DERIV(l*I)=SI<I)*(0.5D0+TI(I)tT) 

10  DERIV(2tI)«TI(I>*(0.5D0+SKI)»S) 

RETURN 
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END 

SUBROUTINE  JACOB ( lELEN » BET JAC » NELEM » NPOI N » COORD » LNODS ) 

IHPLICIT  REAL»8(A-HfO-Z) 

COHMON/UORK/SHAP ( 4 ) » DER I V ( 2 » 4 ) » BNATX ( 3 » 8 ) f POSGP ( 3 ) 
t iWEIGP(3)»OBMAT(3»8)fGPCOD(2»9)fSMATX(3»8t9) 

DIHENSION  COORD< NPOIN 1 2 ) i LNODS ( NELEH » 9) i CARTD( 2 1 4 ) 

REALB8  JAC(2>2)fJACINV(2*2) 

C 

c CALCULATES  THE  JACOBIAN  MATRIX  AND  ITS  DETERMINANT  AND  INVERSE 

C 

DO  5 IDIME-1>2 
DO  S JDIME=1>2 
JAC(IDIME»JDIME)=O.ODO 
DO  5 IN0DE=lf4 
KK=LNODS(IELEM> INODE) 

TEM=COORD<KK»JDIME) 

JAC(IDIME»JDIME)=JAC(IDIME»JDIME)+DERIV(IDIME.INODE)»TEM 
5 CONTINUE 
C 

c CALCULATE  DETERMINANT  AND  INVERSE  OF  JACOBIAN  MATRIX 

C 

DETJAC=JAC(1»1)»JAC(2»2)-JAC(2*1)*JAC(1»2) 

IF(DETJAC.GT.O.O)  GO  TO  10 
HRITE(6»900)  lELEMfDETJAC 

900  FORMAT (//»24HPR0GRAM  HALTED  IN  JACOB  »/fllX» 

»22H  ZERO  OR  NEGATIVE  AREA»/»10Xi 16H  ELEMENT  NUMBER  i I5i 9X»F12.4) 
DO  7 IN0DEsl>4 
MLNODSdELEMflNODE) 

7 WRITE(6>901 )Kf (COORD(Ki J) i J-1 >2) 

901  F0RMAT(/i2X»I3»2F12,4) 

STOP 

C 

10  JACINV(lfl)=JAC(2f2)/DETJAC 
JACINV(2»2)=JAC(lil)/DETJAC 
JAC  I NV  ( 1 » 2 ) »-  JAC  ( 1 » 2 > /DET  JAC 
JACINV(2»1)=-JAC(2»1)/DETJAC 
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C 

c CALCULAT  CARTESIAN  DERIVATIVES  FOR  THE  SHAPE  FUNCTION 

C 

DO  15  IDIHE=lf2 
DO  15  IN0DE=lf4 
CARTD ( IDIHE  f INODE ) =0 . ODO 
DO  15  JDINE=1»2 

CARTD(IDIHE»INODE)=CARTD(IDIME»INODE)+JACINV(IDIMC»JDIME)* 

tDERIV(JDINE>INODE) 

15  CONTINUE 
C 

C CALCULATE  THE  ELEMENT  STRAIN  MATRIX  

C 

N»0 

DO  20  IN0DE=1>4 

M=N+1 

N=M+1 

BMATX ( 1 » M ) =CARTD  < 1 1 INODE ) 

BMATX(liN)=O.DO 
BMATX(2fM)*0.D0 
BM ATX ( 2 » N ) =CARTD ( 2 1 1 NODE ) 

BMATX ( 3 1 M ) =CARTP ( 2 » INODE ) 

BMATX ( 3 » N ) =CARTD ( 1 1 INODE ) 

20  CONTINUE 
RETURN 
END 

SUBROUTINE  STIFPS(NELEM.NPOIN»NGAUSrCOORD»LNODSiBCRT»SP»ROTRfROTI > 

C THIS  ROUTINE  CALCULATES  THE  ELEMENT  STIFFNESS  MATRIX 

IMPLICIT  REALI8(A-H»0-Z) 

REALM  XC(9)f YC<9) f XTf YTf EK(18»18) »BCR(2»18) fYOUNGl >Y0UNG2» 

» P0IS12iG12»THICK*SP(6)»XTIP(2)»YTIP(2)»R0TR(4)»R0TI(4) 

DIMENSION  ESTIF(18>18)fC00RD(NP0INf2)fLN0DS(NELEM>9)>QP(3»3) 
i fBCRT(2>2tl8) 

COMMON/HORK/SHAP ( 4 ) > DER I V ( 2 p 4 ) • BMATX ( 3 > 8 ) f POSGP ( 3 ) 
p pWEICP(3)pQBMAT<3p8)pCPC0D(2p9)pSMATX(3p8p9) 

COMMON/MATRIL/YOUNG1pYOUNG2pPOIS12pG12pIANGLpISTRPCpTHICK 
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COMMON/SUPER/ I SUPR 1 I I SUPR 2 » KE Y 1 f KE Y2 » XT I P . YT I P 

DO  3 I«lf2 
DO  3 J>lf2 
DO  3 K>ifl8 

3 BCRT(IiJfK)sO.ODO 
IF(ISTRPC.E0.2)  GO  TO  5 

C 

C<t«*>CQP3  »8trix  for  islropic  pline  stress  problee....( 
DO  4 I«lf3 
DO  4 J»lf3 

4 OP(I»J)=O.DO 

CONST =Y0UHG1/(1.0-P0IS12*P0IS12)»1.0D0 
QP(l»l)=CONST 
QP(2f2)=CONST 
QP(1»2)=C0NST»P0IS12 
QP(2»1)=CONST»POIS12 
QP(3»3>«CONST*(l.O-P0IS12)/2.O 
GO  TO  6 
C 

C.....C8P3  Matrix  for  anisotropic  plane  stress  problee. 
5 P0IS21=P0IS12tY0UNG2/Y0UNGl 
Q1 1 =Y0UNGl/( 1 . DO-POISl2tPOIS21 ) 
012*Y0UNG2»P0IS12/(1 .DO- P0IS12*P0rS21 ) 
Q22=Y0UNG2/(1.D0-P0IS12»P0IS21) 

066^012 

AHGL = I ANGUDATAN  ( 1 . ODO ) /45 
SIN1=DSIN(ANGL) 

COSl-DCOS(ANGL) 

SIN2»SIN1»SIN1 

C0S2-C0S1»C0S1 

SIN4=SIN2tSIN2 

COS4=COS2>COS2 

SIM3-SIN2l:SINl 

C0S3>C0S2tC0Sl 

QP(1 f 1 )-Ql 1»C0S442« ( Q12t2»Q66 ) *SIN2f C0S2tQ22«SIN4 
QP(l>2)«(QlHQ22-4*Q66)tSIN2*C0S2tQ12*(SIM4IC0S4) 
QP(2>2)=0n>SIN4t2*(012i2tQ66)»SIM2tC052+Q22fC0S4 
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W>(l»3)«(Qll-Q12-2tQ66)»SINltC0S3+(Q12-Q22+2tQ66)»SIN3»C0Sl 

OP(2f3)«(Qll-Q12-2»066)tSIN3»C0Sl+(Q12-Q22+2*Q66)tSINl»C0S3 

QP(3»3)*(01HQ22-2*Q12-2tQ66)»SIN2»C0S2+Q66*(SIN4+C0S4) 

QP(2rl)>QP(if2) 

QP<3fl)>^QP(lf3) 

QP(3.2)=QP(2»3) 

6 CONTINUE 

C LOOP  OVER  EACH  ELEMENT 

C 

DO  70  IELEH= If HELEN 
NDPE=18 

C INITIALIZE  THE  ELEMENT  STIFFNESS  MATRIX 
DO  20  II«lfl8 
DO  20  J>lfl8 
20  ESTIF(IIfJJ)=O.DO 
KGASP=0 

IF(IELEM.E0.ISUPR1.0R.IELEM.EQ.ISUPR2)  GO  TO  68 

C INTER  LOOPS  FOR  AREA  NUMERICAL  INTEGRATION 

DO  50  IGAUS=lfNGAUS 
DO  SO  JGAUS-lfNGAUS 
KGASP=KGASP+1 
EXISP=POSGP(IGAUS) 

ETASP=POSGP(JGAUS> 

C EVALUATE  THE  SHAPE  FUNCTIONS! ELEMENTAL  VOLUMEf  ETC. 

C 

CALL  SHAPE(EXISPfETASP) 

CALL  JACOB ( lELEM  t DET JAC f NELEM f NPOIN  f COORD  f LNODS ) 

DVOLU=  DET JACtUE I GP ( IGAUS ) «UE I GP ( JGAUS ) »THI CK 
C 

C EVALUATE  <QP>LB>  MATRIX 


DO 

10 

II*lf3 

DO 

10 

JJ«lf8 

OBMATCIl 

:f JJ)=0 

DO 

10 

KK=lf3 

QBHAT<IIfJJ)=QBMAT(IIf JJ)-fOP(IIfKK)«BMATX(KKf JJ) 

10  CONTINUE 

C CALCULATE  THE  ELEMENT  STIFFNESSES  (UPPER  TRIANGLE  PORTION) 
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DO  30  II«lf8 
DO  30  JJ>IIf8 
DO  30  K»l>3 

30  ESTIFdlf JJ)«ESTIF(IIf JJ)iBHATX(K»II)«QBHAT(KfJJ)«DVOLU 
50  CONTINUE 

C CONSTRUCT  THE  LOWER  TRIANGLE  OF  THE  STIFFNESS  HATRIX 

DO  60  11=1 f 8 

DO  60  JJ=lf8 

60  ESTIF(JJ>II)=ESTIF(II>JJ) 

60  TO  67 
68  CONTINUE 

IF(IELEM.EQ.ISUPR2)  GO  TO  2A 

KEY=KEY1 

XT=XTIP(1) 

YT=YTIP(1) 

60  TO  23 
24  KEY=KEY2 
XT=XTIP(2) 

YT=YTIP(2) 

23  IF(KEY.EQ.l)  NDPE=IO 
NNPE=NDPE/2 
DO  19  I=lfNNPE 
XC( I ) *COORD( LNODS ( lELEN » I ) » 1 ) 

YC ( I ) =COORD ( LNODS ( I ELEN » I ) f 2 ) 

19  CONTINUE 

URITE(6.40)XTfYTfKEYfNNPE 

40  F0RMAT(/f3X»'XT='»F6.2f7Xf 'YT='fF6.2»7Xi'KEY='»I2»7Xf 'NNPE='fI2»/) 
DO  42  I=1>NNPE 
42  URITE(6>44)I«XC(I)iIfYC(I) 

44  F0RNAT(/f3X.'XC('»Il*')='»F7.3»10X.'YC('»Il»')='»F7.3»//) 
IF(ISTRPC.EQ.l)  GO  TO  55 

CALL  H YBRD2 ( KEY » XC  f YC » XT » YT  f EK  f BCR  f NNPE » SP » ROTR » ROT I ) 

KEY=2 
60  TO  57 

55  CALL  HYBRDl(KEY.XC»YC»XTfYTfEKfBCR»NNPE) 
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K0UNT=2 

IF(IELEH.EO.ISUPRl)  KOUNT=l 
DO  66  K=lfKEY 
DO  66  N-ifNDPE 

66  BCRT(K0UNT>K>N)=6CR(K»N) 

DO  69  I=lrNDPE 

DO  69  Jrl.NDPE 

69  ESTIF(I»J)*EK(IfJ) 

67  «RITE(1)  ESTIF 

70  CONTINUE 
RETURN 
END 

SUBROUTINE  LOADPS ( NELEH  > NPOIN i NGAU5 , COORD  > LNODS  > ELOAD ) 

C =r=rr=====rr===rrrrr=rr===r======r==rrrr=r===r=r=rxxsr= 

C THIS  ROUTINE  IS  TO  CALCULATE  THE  NODAL  FORCES  FOR  EACH  ELEMENT 

IMPLICIT  REAL»8(A-H»0-Z) 

COMMON/HORK/SHAP ( 4 ) i DERI V ( 2 » 4 ) . BMATX ( 3 . 8 ) » POSGP ( 3 ) 
f fMEIGP(3)f0BMAT(3»8)»GPC0D(2»9)»SMATX(3»8»9) 

DIMENSION  PGASH  < 2 ) » DGASH ( 2 ) » NOPRS ( 3 ) r PRE  SS ( 3 f 2 ) t 

» POINT ( 2 ) t STRAN ( 3 ) » STRES ( 3 ) » COORD ( NPOIN . 2 ) » LNODS ( NELEM  f 9 ) i 
• EL0AD(NELEMfl8) 

DO  10  I=1»NELEM 
DO  10  J-l>18 
10  ELOAD(I»J)-O.DO 

C READ  NODAL  POINT  LOADS 

READ(5t930>  NEDGE 
930  FORMATdS) 

HRITE(6f93S)  NEDGE 

935  FORMAT (1 HO »5X. 21 HNO.  OF  LOADING  EDGES=»I5) 

URITE(6>940) 

940  F0RMAT<1H0.5X»38HLIST  OF  LOADED  EDGES  AND  APPLIED  LOADSf/) 
NODEG-2 

C LOOP  OVER  EACH  LOADED  EDGE 

DO  160  IEDCE=lf NEDGE 

C READ  DATA  LOCATING  THE  LOADED  EDGE  AND  APPLIED  LOAD. 

RE AD ( 5 1 945 ) NEASS • ( NOPRS ( I ODEG ) f IODEG-1 > NODEG  > 

945  F0RMAT(3I5) 
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VRITE(6f?50)  NEASSi(N0PRS(I0DEG)iI0DE6=l>N0DEG) 

950  F0RHAT(I10fSXi2I5) 

READ ( 5 f 955 > ( ( PRESS ( I ODEG 1 1 DOFH ) f I ODEG= 1 » NODEG ) 1 1 DOFN= 1 , 2 ) 
HR I TE ( 6 f 956 ) ( ( PRESS ( lODEG  f I DOFN ) » IODEG= 1 f NODES ) » I DOFN= 1 » 2 ) 

955  F0RHAT(4D10.3) 

956  F0RHAT(3X»4F10.3) 

ETASP=-1.D0 

C 

C.... INTER  LOOP  FOR  LINEAR  NUMERICAL  INTEGRATION. 

DO  160  IGAUS=1»NGAUS 
EXISP=POSGP(IGAUS) 

C 

C EVALUATE  THE  SHAPE  FUINCTIONS  AT  THE  SAMPLING  POINTS. 

CALL  SHAPE(EXISP.ETASP) 

C 

C CALCULATE  COMPONANTS  OF  THE  EQUIVALENT  NODAL  LOADS. 

C 

DO  110  ID0FN^lf2 
PGASH(IDOFN)=O.DO 
DGASH(IDOFN)=O.DO 
DO  no  I0DEG=1>N0DEG 

PGASH ( I DOFN ) =PGASH ( I DOFN ) +PRESS  < lODEG  f IDOFN ) tSHAP ( 1 ODEG ) 

1 1 0 DG ASH (I DOFN ) =DGASH (IDOFN) +COORD ( NOPRS ( I ODEG ) f I DOFN ) I 
f DERIV(lflODEG) 

DV0LU=WEI6P(IGAUS) 

PXC0M=^DGASH(1)  tPGASH  ( 2 ) -DGASH  ( 2 ) tPG  ASH  (1 ) 

PYCOM=DGASH ( 1 ) tPGASH ( 1 ) +DG ASH  < 2 ) tPGASH ( 2 ) 

C ASSOCIATE  THE  EQUIVALENT  NODAL  EDGE  LOADS  WITH  AN  ELEMENT,. 

DO  120  IN0DE«lf4 
NLOCA^NODS  (NEASS , INODE ) 

IF(NL0CA.EQ.N0PRS(1)>  GO  TO  130 
120  CONTINUE 
130  JNODE::INODEiNODEG-l 
K0UNT=0 

DO  140  KNODE-INODEiJNODE 

K0UNT=K0UNT+1 

NGASH=(KN0DE-l)t2tl 
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HGASH=(KN0DE-l)t2f2 
IF(KN0DE.GT.4)  NGASH=1 
IF(KN0DE.GT.4)  HGASH=2 

ELOAD (NEASS » NGASH ) =ELOAD ( NE ASS » NGASH ) +SHAP ( KOUNT ) tPXCOM*DVOLU 
1 40  ELOAD ( NEASS . HGASH ) =ELOAD ( NEASS i HGASH ) +SHAP ( KOUNT ) »PYCOM»DVOLU 
160  CONTINUE 
RETURN 
END 

SUBROUTINE  HYBRD 1 < KEY i XC » YC » XT i YT t EK > BCR . NNPE ) 

C cr»=rrr=rsrr==r=rrr====r==rrr===s===rrrrsr==rrr=r=r==r=r= 

COHPLEXtS  ZET»Z»FF2(18)»FF3(18).ZET4»CZZfZETK»ZAfZC»ZEiZDfCI»CZK 
DIMENSION  W(S)>Y<5}>G(18fl8)>XLCN(9)fYLCN(9)»VA(4S)»UB(45)» 
fBCR(2>18)fBK(18fl8)>EKa8fl8)»EKl(171)iUK(63)»Hl(9i9)>H2(9>9)> 
»VINVl(45)»VINV2(45)»TEM(18»18)fVVG(18»18)»XC(9)»YC(9) 
C0HH0N/HATRIL/Y0UNGl>Y0UNG2fP0IS12>G12rIANGLrISTRPC>THICK 
DATA  W/. 23692691 .4786287»0.5688889»0,4786287i0. 2369269/ 

DATA  Y/-0 . 9061799 » -0 . 5384693 »0.0»0.5384693»0. 9061799/ 
CI=CMPLX(O.Ofl.O) 

ETA=(3-P0IS12)/(1+P0IS12) 

SHU=0 . 5* YOUNG 1 / ( 1 +P0 I SI 2 ) 

NT=KEY*9 
DO  1 I=l»45 
VB(I)=0.0 

1 VA(I)=0.0 

NDPE=NNPE»2 
DO  4 I=l»18 
DO  3 K=lt2 

3 BCR(K>I):=0.0 
DO  4 J=liie 
EK(I>J)=0.0 
VVG(I»J)=0.0 

4 G(IfJ)»0.0 

C Rotate  coordinate  to  local  svstea 

ES=SQRT  ( ( XC  ( 1 ) -XT ) »»2+ ( YC  ( 1 ) -YT ) *1:2 ) 

SIN=(YC<1)-YT)/ES 

C0S=(XC(1)-XT)/ES 

XLCT«XT»COS+YT«SIN 
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YLCT*-XTtSIN+YT»COS 
DO  10  IxliNNPE 

XLCN( I )=(XC( I )»C0S+YC( I >*SIN-XLCT)/ES 
YLCM( I ) = { -XC ( I ) »SIN+YC ( I ) »COS- YLCT ) /ES 
10  CONTINUE 

DO  41  ISI=1>4 
UX=YLCN( ISI+1 )-YLCN(ISI ) 

UY=XLCN(ISI)-XLCN(ISI+1) 

DO  41  IIsl«5 

XX=(XLCN(ISI+l)+XLCN(ISI)+(XLCN(ISHl)-XLCN(ISI))»Y(II))/2.0 

YY=(YLCN(ISHl)+YLCN(ISI)+(YLCN(ISI+l)-YLCN(ISI))*Y(II))/2.0 

Z=CMPLX(XX»YY) 

ZET=CSQRT(Z) 

ZET4=Z»Z 

CZZ=C0NJ6(Z) 

ZETK=1.0/ZET4 

KK=1 

IF(KEY.EQ.2)  GO  TO  1011 
DO  1010  K=lf9 
FF2(K)=0.0 
ZETK=ZETK»ZET 
CZK=CONJG(ZETK) 

KK=-KK 

IF(ABS(UY).LE.0.1E-6)G0  TO  1009 

ZD=CZK»CZZ-KKtZETK*Z 

ZC=ZETK»(CZZ-Z) 

FF2 ( K ) = ( 0 . 5*K» ( K-2 ) »ZC+K*ZD ) »UY/2 . 0 

1009  ZE=CZK*CZZ»(CZZ-Z) 

FF3(K)=(ETA»ZETK»ZET4+KK*CZK*CONJG(ZET4)+0.5»K*ZE)/4.0 

IF(ABS(UX).LE.0.1E-6)G0  TO  1010 

IB=K+2+2»KK 

ZA*(K-2)»ZETK-2.0*CZK 

FF2(K)=-K»<CZZ«ZA-IBtZETK»Z)ltCI*UX/4.0+FF2<K) 

1010  CONTINUE 

GO  TO  2000 

1011  DO  1012  Ml>9 
FF2(K)=0.0 
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FF2(K+9)=0,0 

ZETK*ZETK»ZET 

CZK=CONJG(ZETK) 

KK=-KK 

IF(ABS(UY).LE.0.1E-6)C0  TO  1014 

ZD=CZK*CZZ-KK*ZETK*Z 

ZC*ZETK»(CZZ-Z) 

FF2(K)=(0.5m<K-2)tZC+K*ZD)*UY/2.0 
FF2(K+9)=(FF2(K)-K»ZD*UY)»CI 
1014  ZE=CZK»CZZ*(CZZ-Z) 

FF3 ( K ) = ( ETA»ZETK*ZET4+KK»CZK»C0N JG ( ZET4 ) +0 , 5*K*ZE ) /4 . 0 

FF3(K+9)=(FF3<K)-,25»KtZE)»CI 

IF(ABS(UX).LE.0,1E-6)G0  TO  1012 

IB«K+2+2tKK 

ZA=<K-2)«ZETK-2.0tCZK 

FF2(K)»-Kt(CZZ*ZA-IB»ZETK*Z)ICI»UX/4.0+FF2(K) 

IB=K+2-2»KK 

ZA=(K-2)»ZETK+2.0»CZK 

FF2(K+9)=Kt(CZZ*ZA-IB*ZETK*Z)*.25»UX+FF2(K+9) 

1012  CONTINUE 
2000  KJ=0 

DO  41  K=l>9 
L*2»ISI-1 
DO  40  J=lfKEY 
I=J*9+K-9 

G(I»L)  =G(IfL  )+0.5»W<II)*AIMAG(0.5»(1.0-Y(II))»FF2<I)) 
G(I»L+l)=G(IfL+l)+0.5«W(II)»  REAL(0.5»(1,0-Y(II))»FF2(I)) 
G<I.L+2)*G(IiL+2)+0.5*W(II)*AIMAG(0.5*(l,0+Y(II))«FF2(I)) 

40  G(I>L>3)sG(I>L43)IO»Sl:W(II)4;  REAL(0,5I;(1,0+Y<II))*FF2(I)) 

DO  41  J=lfK 

KJ*KJ+1 

VA(KJ)=VA<KJ)+0.5«U(imAIMAG(FF2(K)*FF3(J)+FF2(J)»FF3(K))/SHU 

IF(KEY.EQ.l)  GO  TO  41 

I=K+9 

L=J+9 

VB(KJ)*VB(KJ)+0.5*W(II)*AIMAG(FF2(I)IFF3(L)+FF2(L)*FF3(I))/SMU 

41  CONTINUE 
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IF(KEY.EQ.l)  GO  TO  6A 
DO  500 
6(IIIf2)«0.0 
6(III»l)s2.0*G(IIM) 

II=III+9 

6(II>2)»2.0«6(II>2) 

6(II»1)»0.0 
DO  455  >1>4 
JU=J+1 
JL=10-J 

6(111 »2»JL-l)=G<III»2tJU-l) 

G(IIIf2tJL)=-G(III»2*JU) 

G(II»2»JL-l)=-G(IIf2*JU-l) 

455  6(II»2*JL)»G(II>2«JU) 

500  CONTINUE 

DO  SOI  I:=lrl8 

501  G(11>I)=0»0 

DO  63  I-lf45 
VA(I)=VA(I)*2.0 

63  VB(I)=VB(I)4:2.0 
VB(2)=0.0 
VB(3)=1.0 

DO  62  I«3>9 
62  VB((ItI-I)/2+2)-0,0 

64  CONTINUE 
N»9 

CALL  LINV2P(VA»N»VINVlfIDGT»Dl»D2fWKfIER) 
IF(IER.NE.O)STOP 

IF(KEY.EQ»2)  CALL  LINV2P(VBfNtVINV2»IDGT,Dl»D2fHK»IER) 
IF(IER.NE,0)ST0P 
C Calculate  CHlinv.t  CG3 

C SSSXSXSSSSSZX3SSSS3SS 

DO  111  J-lfNDPE 
DO  no  I»li9 
II=I»(I-l)/2 
BKdf  J)s0.0 
DO  no  K>li9 
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IK=IHK 

IF(K.GT.I)  IK=K*(K-l)/2+I 
no  BK(IrJ)=BK(IrJ)-f0IN01(IK)«G(K>J) 

IF(KEY.EQ.l)  GO  TO  111 
DO  121  I>10>18 

H=I-9 

IJ=Mlt(M-l)/2 
BK(I>J)=0.0 
BO  121  K=lf9 
IK=IJfK 

IF(K.GT.M)  IK=K*(K-l)/2+M 
121  BK(I»J)=BK(If J)+VINV2(IK)tG(K+9»J) 

111  CONTINUE 

C..... Calculate  Eleaent  Stiffness  aatrix  and  Stress  Intensity  Matrix  t 
IJ=0 

DO  112  I=1>NDPE 
DO  112  J=lfl 
IJ=IJ+1 
EK1(IJ)=0.0 
DO  112  K=1»NT 

112  EKl(IJ)=EKl(IJ)+IiK(K»J)»G(K»I)*THICKIYOUNGl 
E=S0RT(2.0/ES) 

DO  113  J=1»NDPE 
BCR(2»J)=0.0 

IF(KEY.EQ.2)  BCR(2f J)=E*BK(10» J) 

113  BCR(1>J)«E«BK(1>J) 

C Rotate  •CBCR>  aatrix  and  eleaent  stiffness  aatrix  to  Global  systeat 
N»0 

DO  3001  I»1>NDPE 
DO  3000  J«1»I 
N=N+1 

TEM(IfJ)=EKl(N) 

3000  TEH(J»I)=EK1(N) 

DO  3001  J=1,KEY 

3001  BK(J>I}sBCR(J>I) 
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DO  3002  >ifKEY 
DO  3002  I-2fNDPEi2 

BCR<J*IMl)*BK<JfIMl)»COS-BK(JfI)tSIN 
3002  BCR(J»I)  *BK<J»IM1)»SIN+BK(J»I)»C0S 
DO  79  J=2fNDPE»2 
I=J-1 

WG(IfI)=COS 
VVG<J»J)=COS 
WG(IfJ)=SIN 
79  VVG(J»I)=-SIN 


DO 

80 

I=1»NDPE 

DO 

80 

J=1»NDPE 

DO 

80 

K=1»NDPE 

DO 

80 

N=1»NDPE 

80  EK(I.J)=EK(If  J)+WG(K.I)*TEM(KiN)»VVG(N.J) 

RETURN 

END 

SUBROUTINE  HYBRD2( KEY » XC » YC » XT » YT i EK » BCR f NNPE » SP » ROTR i ROTI ) 


COMPLEX  ZMUl » ZMU2 » PI » P2  f 01 » 02 » Z1 » Z2 1 ZETl » ZET2 » ZETKl i ZETK2  f 
{ ZTEMl»ZTEM2»ZMUlSiZMU2S 

DIMENSION  W(5)#Y(5)»G(18fl8)fXLCN(9)»YLCN(9)»VA(45)fVB(45)» 
fBCR(2»18)»BKa8»18)»EK(18»18)»EKl(171)»WK(63)»BCRTEM(2»18)f 
»VINVl(45)»VINV2(45)iTEN(18fl8)»VVG(18»18)»XC(9)>YC<9) 
MR0TR(4)»R0TI(4)»SP(6)»FF2(18)»FF3(18)»FF4(18)fFF5(18) 
l>BSTl(18)>BST2(18)>BNDl(18)fBND2(18) 

REAL  Mll>M12fM21fM22 

C0MM0N/NATRIL/Y0UNGl»Y0UNG2fP0IS12»G12fIANGLf ISTRPC»THICK 
DATA  W/ . 2369269 » . 4786287  f 0 . 5688889  r 0 . 4786287 » 0 . 2369269/ 
DATA  Y/-0 . 9061799# -0 . 5384693 . 0 . 0 1 0 . 5384693 » 0 . 9061799/ 
CI=CMPLX(0.0»1.0) 

NT*KEYt9 
DO  1 I=l#45 
VB(I)*0.0 
1 OA(I)sO.O 


NDPE=NNPE«2 


00  4 I=l>18 
DO  3 K=l»2 
BCRTEH(KrI)=0.0 
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3 BCR(K>I)>0.0 
DO  4 >1>18 
BK(IfJ)xO.O 
EK(I>J)=0.0 
VVG(I>J)-0.0 

4 G(I>J):0.0 

CALL  ROOTS(ROTR*ROTI.SP) 

2MU1=CMPLX ( ROTR  < 2 ) > ROT I ( 2 ) ) 

ZMU2=CHPLX ( ROTR  < 4 ) » ROTI ( 4 ) ) 

ZHUiS^ZHUltZHUl 
ZHU2S=ZHU2»ZHU2 
URITE(6f33)ZHUlfZHU2 
P1=SP(1)»ZHU1S+SP(2)-SP(4)»ZHU1 
P2=SP ( 1 ) IZMU2S+SP ( 2 ) -SP ( 4 ) »ZMU2 
Q1=(SP(2)*ZMU1S+SP(3)-SP(5)»ZMU1)/ZMU1 
02= ( SP  < 2 ) 4ZMU2S+SP ( 3 ) -SP ( 5 ) »ZHU2 ) /ZMU2 
33  F0RMAT(//»3X»'ZHUl='»E12.5f ' iSE12,5»10X» 

f 'ZMU2=SE12.5f'  i'»E12.5) 

C 

C Rotate  coordinate  to  local  s»stea 

ES=SQRT ( ( XC ( 1 ) -XT ) *»2+ ( YC ( 1 ) -YT ) »t2 ) 

SINN=(YC(1)-YT)/ES 

C0SS=<XC(1)-XT)/ES 

XLCT=XT*COSS+YT»SINN 

YLCT=-XT»SINN*YT»COSS 

DO  10  I=liNNPE 

XLCN ( I ) » ( XC ( 1 ) 4C0SS+ YC ( I ) »S1 NN- XLCT ) /ES 

YLCN(I)=(-XC(I)»SINN+YC(I)*COSS-YLCT)/ES 
10  CONTINUE 

DO  42  ISI=1>4 

UX=YLCN( ISI+1 )-YLCN( ISI ) 

UY=XLCN(ISI)-XLCN(ISI+1) 

DO  41  II=lfS 
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XX=(XLCN(ISI+l)+XLCN(ISI)+(XLCN(ISI+l)-XLCN(ISI))*Y(II))/2.0 
YY=(yLCN(ISI+l)+YLCN(ISI)+(YLCN(ISI+l)-YLCN(ISI))»Y(II))/2.0 
THETAl=ATAN2(R0TK2)m  » (XX+R0TR(2)»YY) ) 
THETA2*ATAN2(R0TI(4)m  » (XX+R0TR(4)m) ) 

R1 =SQRT ( ( XX+ROTR ( 2 ) *YY ) »*2+ ( ROT I ( 2 ) » Y Y ) *42 ) 

R2=SQRT( ( XX+ROTR ( 4 ) »YY ) *42+ ( ROT I ( 4 ) 4YY ) 442 ) 
XR1=R14C0S(THETA1) 

YI1=R14SIN(THETA1) 

XR2=R24COS(THETA2) 

YI2=R24SIN(THETA2) 

Z1=CMPLX<XR1»YI1) 

Z2=CHPLX(XR2iYI2) 

ZET1=CSQRT(Z1> 

ZET2=CSQRT(Z2) 

ZETK1=1.0/Z1 

ZETK2=1.0/Z2 

KK=1 

IF(KEY.EQ.2)  GO  TO  1011 

DO  1010  K=l»9 

FF2(K)=0.0 

FF4(K)=0.0 

ZETK1=ZETK14ZET1 

ZETK2=ZETK24ZET2 

ZTEM1*ZETK14Z1 

ZTEH2=ZETK24Z2 

KK=-KK 

IF(KK.EQ.l)  GO  TO  13 
H11=-R0TI(2)/R0TI(4) 

M12= ( ROTR ( 4 ) "ROTR ( 2 ) )/ROTI ( 4 ) 

M21=0.0 
H22=-1.0 
60  TO  14 
13  Hlls-1.0 
H12=0.0 

«21=(R0TR(2)-R0TR(4))/R0TI(4) 

M22=-R0TI(2)/R0TI(4) 


14  CONTINUE 
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C J*ROTR (2 ) »REAL ( ZETKl ) + ( ROTR ( A ) *M1 1-ROT I ( 4 ) »H21 ) *REAL ( ZETK2 ) 

I -R0TI(2)*AIHAG(ZETK1)-(R0TR(4)»M2HR0TI(4)»H11)*AIMAG(ZETK2) 
IF(ABS(UY).LE.0.1E-A)G0  TO  1009 
B J=REAL  ( ZETK  mHl  1 *RE  AL  ( ZETK2 ) - M2 1 » A I HAG  ( ZETK2 ) 

FF2(K)»-KtCJ»UY 
FF4<K)=  KtBJ»UY 

1009  FF3(K)-REAL(P1)|[REAL(ZTEH1)-AIHAG<P1)»AIHAG(ZTEH1) 

t 4Hllt(REAL(P2)*REAL(ZTEH2)-AIMAG(P2)*AIMAG(ZTEH2) ) 

I -M21»(REAL(P2)*AIMAG(ZTEM2)+AIMAG(P2)*REAL(ZTEM2) ) 

FF5<K)=REAL(Q1)»REAL(ZTEH1)-AIMAG(Q1)*AIMAG(ZTEH1) 

« 4H11«(REAL(Q2)«REAL(ZTEH2)-AIHAG(Q2)»AIHAG(ZTEH2) ) 

♦ “M21t<REAL(Q2)»AIHAG(ZTEM2)+REAL(ZTEH2)*AIHAG(Q2) ) 

IF(ABS(UX).LE.0.1E-6}G0  TO  1010 

AJ=REAL(ZMU1S)»REAL<ZETK1)-AIMAG(ZHU1S)*AIMAG(ZETK1) 

I +H11»(REAL(ZMU2S)*REAL(ZETK2)-AIMAG(ZMU2S)*AIMAG(ZETK2)) 

% “M21»(REAL(ZMU2S)»AIHAG(ZETK2)+REAL(ZETK2)tAIHAG(ZMU2S) ) 
FF2(K)=K*AJtUX+FF2(K) 

FF4{K)=-K»CJ»UX+FF4(K) 

1010  CONTINUE 

60  TO  2000 

1011  DO  1012  K>1>9 
FF2(K)=0.0 
FF2(K+9)=0,0 
FF4(K)=0.0 
FF4(K+9)=0.0 
ZETK1=ZETK1»ZET1 
ZETK2=ZETK2»ZET2 
ZTEMUZETKUZl 
ZTEM2=ZETK2»Z2 
KK*-KK 

IF(KK.EQ.l)  GO  TO  15 
M11*-R0TI(2)/R0TI(4) 

M12»(R0TR(4)-R0TR<2))/R0TI(4) 

H21«0.0 
M22»-1.0 
60  TO  16 
15  Mll*-1.0 


H12>0.0 

H21=(R0TR<2)-R0TR(4))/R0TI(4) 

M22-R0TI(2)/R0TI(4) 
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16  CONTINUE 

CJ=R0TR(2)tREAL(ZETKl)+<R0TR<4)*Mll-R0TI(4)»N21)*REAL(ZETK2) 

• -ROT I ( 2 ) lAIMAG ( ZETK 1 ) - ( ROTR ( 4 ) IM2 1 +ROT I ( 4 ) *M1 1 ) »A IMAG ( ZETK2 ) 
C JN= ( ROTR ( 4 ) »H1 2-ROT I ( 4 ) »M22 ) «RE AL ( ZETK2 ) -ROT I ( 2 ) »RE AL ( ZETKl ) 

I - < ROTI ( 4 ) »N1 2+ROTR ( 4 ) *M22 ) 4AIHAG ( ZETK2 ) 

% -R0TR(2)tAIHAG(ZETKl) 

IF(ABS(UY).LE.0.1E-6)C0  TO  1014 

BOREAL  ( ZETKl ) 4H1 1 CREAL  ( ZETK2 ) -H21  »AIHAG  ( ZETK2 ) 

B JN=H12*REAL ( ZETK2 ) -AIHAG ( ZETKl ) -H22»AINAG( ZETK2 ) 

FF2(K)  «-K»CJ  *UY 
FF2(K+9)=-K$CJN»UY 
FF4(K)  s K4BJ  »UY 
FF4(K+9)=  KtBJNtUY 

1014  FF3(K)  -REAL(P1)«REAL(ZTEN1)-AINAG(P1)«AINAG(ZTEH1) 
t f HI 1»( REAL (P2 ) >REAL ( ZTEH2 ) -AIHAG ( P2 ) tAIHAG ( ZTEH2 ) ) 

I -H2 1 » ( REAL ( P2 ) *A I HAG ( ZTEH2 ) ♦ A IHAG ( P2 ) »RE AL ( ZTEH2 ) ) 

FFS(K)  =REAL<Ql)tREAL(ZTEHl)-AIHAG(Ql)»AIHAG(ZTEMl) 

% +Hll»(REAL(Q2)l:REAL(ZTEH2)-AIMAG(Q2)*AIHAG<ZTEH2) ) 
i -H21 I ( REAL (02 ) CAIHAG ( ZTEH2 ) +RE AL ( ZTEM2 ) CAIHAG ( 02 ) ) 
FF3(K+9)=H12»<REAL<P2)*REAL<ZTEH2)-AIHA3<P2)»AIMAG(ZTEH2)) 
i -H22»(REAL<P2)*AIHAG(ZTEH2)+REAL(ZTEH2)*AIHAG(P2) ) 

» -REAL(Pl)tAIHA6(ZTEHl)-REAL(ZTEHl)*AIHA6(Pl) 

FF5 ( K49 ) =H1 2t ( REAL ( 02 ) «REAL ( ZTEH2) -AIHAG ( 02 ) 4AIHAG  < ZTEH2 ) ) 

I -H22«(REAL(02)IAIHA6(ZTEH2)4REAL(ZTEH2)«AIHA6(02)) 

t -REAL(01)*AIHAG(ZTEH1)-REAL(ZTEH1)»AIHAG(Q1) 

IF(ABS(UX).LE.0.1E-6)G0  TO  1012 
AJ=REAL(ZHUlS)tREAL(ZETKl)-AIHAG(ZHUlS)*AIHAG(ZETKl) 

• 4H1 1 1 ( REAL ( ZHU2S  > 4REAL ( ZETK2 ) -AIHAG ( ZHU2S ) »AI HAG ( ZETK2 ) ) 

% -H21»(REAL ( ZHU2S) 4AIHAG ( ZETK2 ) +REAL  < ZETK2 ) lAIHAG ( ZHU2S) ) 

AJN=H12» ( REAL ( ZHU2S ) »REAL ( ZETK2)-AIHAG ( ZHU2S) tAIHAG ( ZETK2) ) 

• -H22» ( REAL ( ZHU2S ) »AI HAG ( ZETK2 ) +RE AL  < ZETK2 ) 4AIHAG ( ZHU2S ) ) 

I -REAL<ZHU1S>*AIHAG(ZETK1)-REAL(ZETK1)»AIHAG(ZHU1S) 

FF2(K)  « K»AJ  *UX+FF2(K) 

FF2(K+9)»  K*AJN»UX+FF2(K+9) 
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FF4(K)  =-K»CJ  *UX+FF4(K) 

FF4(K+9)=-K»CJN*UX+FF4(K+9) 

1012  CONTINUE 
2000  KJ=0 

DO  41  K»lf9 
L*2»ISI-1 
DO  40  J=liKEY 
I=J*9+K-9 

G(IfL)  *G(IfL  )40.5»W(II)»(0.5*(1.0-Y(II))»FF2(I)) 
G(I»L+l)=G(IfL+l)+0.5*W(II)*(0.5»(1.0-Y(II))»FF4(I)) 
G(IiL+2)»G(I»L+2)+0.5«W(II)*(0.5»(1.04Y(II))*FF2(I)) 

40  G(IfL+3)=G(I»L+3)+0.5»W(II)»(0,5»(1.0+Y(II))»FF4(I)) 
C 

DO  41  J^lfK 
K>KJ+1 

VA(KJ)=VA(KJ)40.5tW(II)*(FF2(K)»FF3(J)+FF2(J)»FF3(K) 
t +FF4(K)*FF5(J)+FF4(J)*FF5(K)) 

IF(KEY.EQ.l)  GO  TO  41 

I«K+9 

L=J+9 

VB(KJ)=VB(KJ)+0.5*W(II)*(FF2(I)IFF3(L)+FF2(L)*FF3(I) 
« +FF4(I)tFF5(L)+FF4(L)»FF5(D) 

41  CONTINUE 

42  CONTINUE 
IF(KEY.EQ.l)  GO  TO  44 
DO  500  III=li9 
G(III>2)=0.0 
G(IIM)-2.0»G(III>1) 

II=III+9 

G(II>2)«2.0>G(II»2) 

G(IIfl)=0.0 
DO  455  >1>4 
JU=J+1 
JL=10-J 

G(III»2tJL-l)=G(IIIf2«JU-l) 

G(III>2tJL)=-G(III>2>JU) 

G ( 1 1 f 2» JL-1 )»-G( 1 1 . 2» JU-1 ) 
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455  6(II>2tJL)>G(II»2tJU) 

500  CONTINUE 

DO  501  1^1 >18 

501  G(11>I)»0.0 

DO  63  I>1>45 
VA(I><:0A(I>t2.0 

63  VB(I)=0B(I)«2.0 
VB<2)=0.0 
VB(3)»1.0 

DO  62  I»3>9 
62  VB((I»I-I)/2+2)=0.0 

64  CONTINUE 
N=9 

CALL  LINV2P(VA>N>VINUl>IDGT>Dl>D2>UKfIER) 
IF(IER.NE.0)WRITE(6>82) 

IF(IER.NE.0)ST0P 

82  FORMAT (//>9X> 'CHECK  MATRIX  Hl'>//) 

IF(KEY.EQ.2)  CALL  LIN02P(UB>Nf VIN02> IDGTrDl >D2>UK> lER) 
IF(IER.NE.0)URITE(6>83) 

83  FORMAT (//>9X I 'CHECK  MATRIX  H2'>//) 

IF(IER.NE.O)STOP 

C Calculate  CHIinv.t  CG3 

C ===================== 

DO  111  >1>NDPE 
DO  no  I»l>9 
II=It(I-l)/2 
BK(IfJ)=0.0 
DO  no  K>1>9 
IK=II+K 

IF(K.BT.I)  IK=K»(K-l)/2+I 
no  BK(I>J)-BK(IfJ)iOINVKIK)ltG(K>J) 

IF(KEY.EQ.l)  GO  TO  111 
DO  121  I>10>18 

M=I-9 

IJ=Mt(M-l)/2 
BK(I>J)=0.0 
DO  121  K»l>9 


IF(K,OT,M)  IK=K*(K-l)/2+H 
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121  BK(IfJ)=BK(I»J)+VINV2(IK)*G(K+9»J) 

111  CONTINUE 

C Cilculat*  Eleaent  Stiffness  aatrix  and  Stress  Intensity  Matrix  1 

IJ=0 

DO  112  I«1>NDPE 
DO  112  J>1»I 
I>IJ+1 
EKKIJ)-0.0 
DO  112  K^lfNT 

112  EKl(IJ)=EKl(IJ)fBK(K>J)tG(KfI)«THICK 
E*SQRT(2.0/ES) 

T1«1-R0TI(2)/R0TI(4) 

T2= ( ROTR ( 4 ) -ROTR ( 2 ) ) /ROT I ( 4 ) 

T3=R0TR (2 ) -ROTR( 4 ) *R0TI ( 2) /ROTI ( 4 ) 

T4=R0TI(4)-R0TI(2)+R0TR(4)*T2 
DO  113  J=lfNDPE 
BSTl(J)sEtTl»BK(l>J) 

BND1(J)=-E»T3»BK(1.J) 

IF(KEY.EQ.l)  GO  TO  23 
BST2(J)»  E*T24BK(10fJ) 

BND2(J)=-EI!T44:BK(10.J) 

BCR(1.J)=BST1(J>+BST2(J) 

BCR(2fJ)>BNBl(J)4BN02(J) 

60  TO  113 

23  BCR(1>J)=BST1(J) 

BCR(2rJ>=BNDl(J) 

113  CONTINUE 

C Rotate  <BCR>  aatrix  and  eleaent  stiffness  aatrix  to  Global  systea! 
N«0 

DO  3001  I»1>NDPE 
DO  3000  J^lrl 
N=N41 

TEH(I>J)=:EK1(N> 
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3000  TEM(JiI)»EKl(N) 

DO  3001  J=l>2 

3001  BCRTEH(J»I)=BCR(J*I) 

DO  3002  J-l>2 

DO  3002  Is2fNDPEf2 

BCR(JiIHl)=BCRTEM(J»IHl)*COSS-BCRTEM(J»I)»SINN 

3002  BCR(J»I)  =BCRTEM<J»IMl)»SINN+BCRTEM(JfI)*COSS 
DO  79  J=2fNDPE»2 

I=J-1 

W6(IiI)=C0SS 

WG(J»J)=COSS 

VV6(I»J)=SINM 

79  VVG(J»I)=-SINN 
DO  80  I=1»NDPE 
DO  80  J=1>NDPE 
DO  80  K=lfNDPE 
DO  80  N=1>NDPE 

80  EK(IfJ)=EK(I»J)+VVG(K»I)*TEM(K»N)*WG(NfJ) 

RETURN 

END 

SUBROUTINE  ROOTS (RORfROIfSPP) 

DIMENSION  XC0(5)»SPP(6)»C0(5)fR000TR(4)fR000TK4)fR0R(4)»R0I(4) 
C0MM0N/MATRIL/Y0UNGl»Y0UNG2fP0IS12fG12»IANGL»ISTRPCf THICK 
H=4 

S11*1/Y0UNG1 

S12*-P0ISl2tSll 

S22*1/Y0UN62 

S66»l/G12 

C Calculate  the  coefficients  of  the  charactaristic  eouation. 

IF(IANGL.EQ.0.0R.IANGL.EQ.90)  CO  TO  2 
ANGL« I ANGLtAT AN (1 . 0 ) /45 
SIN1«SIN(AN6L) 

C0S1»C0S(AN6L) 

SIN2=SINi»SINl 

C0S2»C0S1*C0S1 

SIN4sSIN2»SIN2 
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COS4=COS2*COS2 

SIN3=SIN2*SIN1 

C0S3>C0S2»C0S1 

SPP(l)«SlltC0S4+(2*S12+S66)*SIN2»C0S2+S22tSIM4 
SPP(2)=S12»(SIN4+C0S4)+(S1HS22-S66)»SIN2»C0S2  - 
SPP (3)»S1 USIN4+ ( 24S12+S66 ) *SIN2»C0S2+S22»C0S4 

SPP(4)=(2tSll-2tS12-S66)»SINUC0S3-(2tS22-2*S12-S66)»SIN3*C0Sl 

SPP(5)=(2*S11-2»S12-S6&)*SIN3*C0S1-(2»S22-2*S12-S66)*SIN1*C0S3 

SPP(6)=2t(2»SU+2»S22-4»S12-S66)»SIN2»C0S2+S66*(SIN4+C0S4) 

yRITE(6f8)(SPP(I>fI=l»6) 

8 F0RMAT(2Xf'SPP(I)='»6E14.6»//) 

XC0(1)=SPP<3) 

XC0(2)=-2*SPP(5) 

XC0<3)=2#SPP(2)+SPP(6) 

XC0(4)=-2tSPP(4) 

XC0(5)=SPP(1) 

60  TO  6 
2 XC0(2)=0.0 

XC0(3)*2»S12+S66 

XC0(4)=0.0 

SPP(2)*S12 

SPP(4)=0.0 

SPP(5)=0.0 

SPP(6)=S66 

IF(IAN6L.E0.90)  GO  TO  4 

sppa)=sn 

SPP(3)=S22 
XC0(i)=S22 
XC0(5)=S11 
60  TO  6 
4 XCO<i)sSll 
XC0(5)sS22 
SPP(1)=S22 
SPP(3)=Sll 

6 WRITE<6fl2)(XC0(I)fI-lf5) 

12  F0RMAT(5Xf'XC0(I)='»5E14.6»//) 

Ctt»««Calculate  iht  roots  of  the  charactaristic  eouation# 
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CALL  POLRT ( XCO»  CO » M » ROOOTR » ROOOTI » lER ) 

DO  7 I>1>4 

7 yRITE(6>10)  I>ROOOTR(I)fIfROOOTI(I) 

10  F0RMAT(5Xf 'ROOTR('fIlf')s'tF10.5f6Xi'ROOTI('»Ilf 
DO  14  I«li4 
ROI(I)eROOOTKI) 

IF(IAN6L.EQ.0.0R.IANGL.EQ»90)  60  TO  13 
RORa)=ROOOTR(I) 

CO  TO  14 

13  R0R(I)=0.0 

14  CONTINUE 
RETURN 
END 

SUBROUTINE  POLRT ( XCOF » COF i M f ROOTR i ROOT I » lER ) 

DIMENSION  XCOF ( 1 ) t C0F(1 ) »ROOTR( 1 ) iROOTI ( 1 ) 

IFIT=0 

N=M 

IER=0 

IF<XC0F<N+l))10»25fl0 
10  IF(N)  lSflS*32 
IS  lER^l 
20  RETURN 
25  IER=4 
GO  TO  20 
30  IER=2 
GO  TO  20 

32  IF(N-36)  35»35>30 
35  NX«N 
NXX=N*1 
N2=l 

KJl  = N+1 
DO  40  L^lfKJl 
HT»KJ1-L+1 
40  COF(MT)=XCOF(L) 

45  X0=» 005001 01 
Y0=0. 01000101 
IN=0 
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50  X=X0 

X0=-10,0*Y0 
Y0=-10.0»X 
X=X0 
Y«Y0 
IN=IN+1 
GO  TO  59 
55  IFIT=1 
XPR=X 
YPR=Y 

59  ICT=0 

60  UX=0.0 
UY=0.0 
V =0.0 
YT=0.0 
XT=1.0 
U=C0F(N+1) 

IF(U)  65fl30>65 

65  DO  70  1=1 »N 
L =N-I+1 
TEMP=COF(L) 
XT2=X»XT-Y*YT 
YT2=X»YT+Y*XT 
U=U+TEMP»XT2 
V=V+TEHP»YT2 
FI=I 

UX=UX+FI»XmEMP 
UY=UY-FI»YT*TEMP 
XT=XT2 
70  YT=YT2 

SUMSQ=UX»UX+UY»UY 
IF(SUMSQ)  75fll0f75 
75  DX=(V»UY-U»UX)/SUMSQ 
X«X+DX 

DY»- ( U»U Y+V tux ) /SUMSQ 
Y»YfDY 


78  IF<ABS(DY)tABS(DX)-1.0E-0S)  lOOfSOiSO 
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80  ICT»ICT+1 

IF<ICT-500)  60»85f85 
85  IF(IFIT)100i90»100 
90  IF(IN-5)  50f95f95 
95  IER=3 
GO  TO  20 

iOO  DO  105  L=1>NXX 
HT=KJ1-L+1 
TEMP=XCOF(HT) 

XCOF(MT)=COF(L) 

105  COF(L)=TE«P 
ITE«P=N 
N=NX 

NX=ITEMP 

IF(IFIT)  120i55»120 
no  IF<IFIT)  115»50.115 
115  X=XPR 
Y=YPR 
120  IFIT=0 

122  IF(ABS(Y)-1.0E-4*ABS(X))  135fl25»125 
125  ALPHA=X+X 
SUMSQ=X*X+Y»Y 
N=N-2 
GO  TO  140 
130  X=0.0 
MX=NX-1 
NXX=NXX-1 
135  Y=0.0 
SUH50=0.0 
ALPHA=X 
H=N-1 

140  C0F(2)=C0F(2)+ALPHA*C0F(1) 

145  DO  150  L-2>N 

150  C0F(L+1)=C0F(L+1)+ALPHA«C0F(L)-SU«SQ*C0F(L-1) 
155  R00TI(N2)=Y 
R00TR(N2)»X 
N2-N2F1 
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IF(SUHSQ)  160fl65>160 
160  Y=-Y 

SUHSQ>0.0 
GO  TO  1S5 
165  IF(N)  20f20»65 
END 

SUBROUTINE  FRONT <NELEM»NPOIN»NVFIX»NTOTV»NGAUS»COORDiLNODSf 
f ELOADi ASDIS » NOFIX » IFPRE  »PRESC  t BCRT » MFRON » CRD » SIF ) 

IMPLICIT  REAL»8(A-HiO-Z) 

DIMENSION  C00RD(NP0INf2)»LN0DS(NELEMf?)fEL0AD(NCLEMfl8)» 
fQV(18)»ASDIS(NT0TV)»N0FIX(NVFIX)»IFPRE(NVFIX>2).PRESC(NVFIX»2) 
« fCRD(NP0INf2)f0Ul(18)>SIF(4) 

REALM  XTIP(2)fYTIP(2) 

COMHON/WORK/SHAP ( 4 ) » DERI V ( 2 * 4 ) » BMATX ( 3 » 8 ) i POSGP ( 3 ) 
t ,HEICP(3)»QBMAT(3»8)»GPC0D(2»9)»SMATX(3i8»9) 

COMMON/SUPER/ 1 SUPRl » I SUPR2 » KE Y1 1 KE Y2  f XT I P . YT I P 
C 

DIMENSION  F I XED ( 368 ) » EQUAT ( 60 ) f VECRV ( 60 ) » GLOAD ( 60 ) » LOCEL ( 1 8 ) » 
»GSTIF( 1830) f ESTIF( 18»18> » IFFIX(368) »NACVA(60) ,NDEST( 18) » 

» BCRT(2>2>18) 

NFUNC(I#J)=IF(J»(J-l))/2 

C 

NSTIF=MFR0Nt(MFR0N+l)/2 

C 

C Inirpret  fixity  data  in  vector  

C 

DO  100  I=1»NT0TV 
IFFIX(I)«0 
100  FIXED(I)=O.DO 

DO  no  I»1»NVFIX 
NL0CA=2»(N0FIX(I)-1) 

DO  no  J=l»2 
N6ASH=NL0CAFJ 
IFFIX(NGASH)=IFPRE(I»J) 
no  FIXED(NGASH)=PRESC(I>J) 

C 
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C«»»t»Ch8n4e  the  sifln  of  last  aPF'earsnee  of  each  node 
DO  140  I>1»NP0IN 
KLAST-0 

00  120  lELEh-liNELEH 
NN0DE=4 

IF( lELEH.EQ. ISUPRl . AND.KEYl .EQ. 1 )NN0DE-5 
IF ( lELEH . EQ » ISUPRl . AND . KE Y 1 ,EQ . 2 ) NN0DE=9 
IF ( lELEH . EQ . ISUPR2 ) NN0DE=9 
DO  120  IN0DE=liNN0DE 
IF(LNODS(IELEM»INODE).NE.I)  GO  TO  120 
KLAST=IELEH 
NLAST^INODE 
120  CONTINUE 

IF ( KLAST . NE . 0) LNODS ( KLAST » NLAST ) =- 1 
140  CONTINUE 
C 

C I NITIALIZATIO  N 

DO  150  ISTIF=1»HSTIF 
150  GSTIF<ISTIF)=O.DO 

DO  160  IFR0N=1»NFR0N 
GLOAD(IFRON)=O.DO 
EQUAT(IFRON)=0.00 
VECRV(IFR0N)=0,D0 
160  NACVA(IFR0N)=0 
C 

C PREPARE  FOR  DISC  READING  AND  WRITING  OPERATIONS 

REWIND  1 
REWIND  2 
C 

C ENTER  MAIN  ELEMENT  ASSEMBLY  - REDUCTION  LOOP,.,. 

NFR0N=0 

KELVA-0 

DO  380  I»1>NELEM 
NN0DC«4 

IF( I. EQ. ISUPRl. AND. KEY1,EQ.1)NN0DE=5 
IF ( I .EQ. ISUPRl .AND.KEYl .EQ. 2) NN0DE=9 
IF ( I . EQ . ISUPR2 ) NNODEs? 
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NEVAB=NN0DEt2 
KEVAB==0 
READ(l)  ESTIF 
DO  165  IEVAB=1»NEVAB 
165  LOCEL<IEVAB)»0 

DO  170  J-1>NN0DE 
DO  170  Ksl>2 
NP0SI=(J-1)I2+K 
LOCNO-LNODS(IfJ) 

IF ( LOCNO .GT .0 ) LOCEL ( NPOSI ) = ( LOCNO  -1 ) *2+K 
IF(LOCNO.LT.O)  LOCEUNPOSI ) = (L0CN0+1  )»2-K 
170  CONTINUE 
C 

C START  BY  LOOKING  FOR  EXISTING  DESTINATIONS 

C 

DO  210  IEVAB=lfNEUAB 
NIKNO=IABS(LOCEL(IEOAB) ) 

IF(NIKNO.EQ.O)  GO  TO  210 
KEXIS^O 

IF(NFRON.EQ.O)  60  TO  181 

DO  180  IFR0N=1*NFR0N 

IF(NIKNO.NE.NACVA(IFRON))  GO  TO  180 

KEVAB=KEVAB+1 

KEXIS-1 

NDEST(KEVAB)=IFRON 

180  CONTINUE 

181  CONTINUE 
IF(KEXIS.NE.O)  GO  TO  210 

C 

C NE  NOW  SEEK  NEW  EHPTY  PLACES  FOR  DESTINATION  VECTOR 

C 

DO  190  IFRON=liHFRON 

IF<NACVA(IFRON).NE.O)  60  TO  190 

NACVA(IFRON)=NIKNO 

KEVAB>KEVABtl 

NDEST(KEUAB)=IFRON 

GO  TO  200 
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i?0  CONTINUE 
C 

C THE  NEU  PLACES  NAY  DEHAND  AN  INCREASE  IN  CURRENT  FRONTHIDTH 

C 

200  IF(NDEST(KE0AB).6T.NFR0N)  NFRON=NDEST(KEVAB) 

210  CONTINUE 

C««*»«Asse»ble  element  loads 

DO  240  II^liNEUAB 
IFaOCEL(II).EQ.O)  GO  TO  240 
IDEST=NDEST(II) 

GLOAD ( I DEST ) *GLOAD  < I DEST ) +EL  OAD ( I » 1 1 ) 

C Asseable  the  eleeent  stiffnesses  but  not  in  resoIution» •••#>« > 

DO  220  JJslfll 
JDEST=NDEST(JJ) 

NGASH=NFUNC( IDEST  t JDEST ) 

NG I SH=NFUNC ( JDEST i I DEST ) 

IF ( JDEST , GE . I DEST ) GST I F ( NG ASH ) =GST I F ( NGASH ) +EST I F ( 1 1 » J J ) 

IF( JDEST.lt. IDEST)GSTIF(NGISH)=GSTIF(NGISH)iESTIF(II»JJ) 

220  CONTINUE 
240  CONTINUE 

C. . . . .Re-E>:aein  each  eleeent  node  to  enauire  which  can  be  elieinated. . . 
DO  370  II-l.NEUAB 
NIKNO=  -LOCEL<II) 

IF(NIKNO.LE.O)  CO  TO  370 

C Find  positions  of  variables  readw  for  elieination. .......... . 

DO  350  IFRON^lfNFRON 
IF(NACVA(IFRON).NE.NIKNO)GO  TO  350 

C Extract  the  cofficients  of  the  new  eouation  for  elieination. . . 

DO  250  JFRON=liHFRON 

IF( IFRON .LT . JFRON ) NLOCA=NFUNC ( IFRON . JFRON ) 

IF ( IFRON . CE . JFRON)  NLOCA^NFUNC ( JFRON i IFRON ) 

EQUAT  < JFRON ) =GSTI F ( NLOCA ) 

250  GSTIF(NLOCA)>O.DO 

C.....And  extract  the  corresponding  ridht  hand  sides 

EQRHS-GLOAD( IFRON) 

6L0AD( IFRON )«0. DO 
KELUA«KELVAfl 


186 


Write  eouetions  to  DISC  or  to  TAPE 
WRITE(2)  EQUATf  EQRHSf  IFRON*  NIKNO 

C»««»<De3l  with  eivot 

PIVOT=EQUAT<IFRON) 

EQUAT(IFRON)sO.DO 

C Enauire  whether  present  verieble  is  free  or  prescribed 

IFdFFIX(NIKNO).EQ.O)  60  TO  300 

C.....De8l  with  s prescribed  displacement..... 

DO  290  JFR0N=1>NFR0N 

290  GLOAD ( JFRON) =GL0AD ( JFRON ) -FIXED ( NIKNO ) «E0UAT  < JFRON ) 

60  TO  340 

C Eliminate  a free  variable  - deal  with  the  riaht  hand  side  first... 

300  DO  330  JFR0N=liNFR0N 

GLOAD ( JFRON ) =GL0AD  < JFRON ) -EQUAT ( JFRON ) »EQRHS/PI VOT 

C Now  deal  with  the  coefficients  in  core..... 

IF(EQUAT(JFR0N).EQ.0.0)  GO  TO  330 
NLOCA-NFUNC( Of JFRON) 

DO  310  LFR0N=1> JFRON 

NGASH-LFRONfNLOCA 

310  GSTIF ( NGASH ) =GST I F (NGASH ) -EQUAT ( JFRON ) »EQUAT ( LFRON ) /PIVOT 

330  CONTINUE 

340  EQUAT (IFR0N)=PIV0T 

C..... Record  the  new  vacant  space*  and  reduce  frontwidth  if  possible..... 
NACVA(IFR0N)»0 
60  TO  360 

C. ... .Complete  the  element  loop  in  the  forward  elimination...... 

350  CONTINUE 

360  IF(NACVA(NFR0N).NE.0)60  TO  370 
NFR0N=NFR0N-1 
IF(NFRON.GT.O)  GO  TO  360 
370  CONTINUE 
380  CONTINUE 

C Enter  back-substitution  phase*  loop  backwords  throuSth  variables... 

DO  410  IELVA-1*KELVA 

C.....Read  a new  eouation 

BACKSPACE  2 

READ(2)  EQUAT*  EQRHS*  IFRON*  NIKNO 
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BACKSPACE  2 

C»»*««PreP3re  to  back  - substitute  froo  the  current  esuetion. « • • • « • 
PIVOT«EQUAT(IFRON) 

IF(IFFIX(NIKNO).EQ.l)  VECRV(IFRON)=FIXED(NIKNO) 
IF(IFFIX(NIKNO)»EQ.O)  EQUAT(IFR0N)=0.D0 

Ctot^Bsck-substitute  in  the  current  eou3tioni.<««>> 

DO  400  JFR0N=lfHFR0N 
400  EQRHS=EQRHS-VECRO  < JFRON ) tEOUAT ( JFRON ) 

C*»«>»Put  the  final  values  where  thew  belons. 

IF ( IFFIX ( NIKNO ) . EQ .0 ) VECRV( IFRON ) =EQRHS/PI VOT 
IF(IFFIX<NIKN0).EQ.1)  FIXED(NIKN0)=-E0RHS 
ASDI S ( NIKNO ) =VECRV ( IFRON ) 

410  CONTINUE 

WRITE(6fP00) 

900  F0RHAT(//fl5X»13HDISPLACEMENTS) 

MRITE(6f905) 

905  FORMAT ( 1 HO  * 5X  f 4HN0DE 1 6X » 7HX-D I SP . » 7X » 7HY-DI SP . ) 

DO  450  IP0IN=1»NP0IN 

NGASH=IP0IN«2 

NGISH=N6ASH-1 

CRD(IPOINfi)=COORD(IPOIN>l)-ASDIS(NGISH) 

CRD(IP0IN>2)=C00RD(IP0INf2)-ASDIS(NGASH) 

450  URITE(6f920)IP0INt<ASDIS(IGASH)fIGASH=NGISHfNGASH)>(CRD(IP0INf 
t M)iM=l»2> 

920  F0RHAT(I10f2E14.6f2F15.6) 

URITE(6>925) 

925  F0RMAT(lH0fl5Xf9HREACT10NS) 

HRITE(6f935) 

935  FORMAT ( IHO » 5X » 4HN0DE » 5X  * 7HX-F0RCE » 7X » 7HY-F0RCE ) 

DO  510  IPOIN-lfNPOIN 
W-0CA=(IP0IN-1)*2 
DO  490  ID0FN«lf2 
NGUSH>NLOCAHDOFN 
IF(IFFIX(N6USH).GT.O)  CO  TO  500 
490  CONTINUE 
60  TO  510 

500  NCASH=NL0CAt2 
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NGISH=NLOCAH 

WRITE(6i945)  IPOINi  (FIXED(IGASH)»IGASH=NGISH»NBASH) 

510  CONTINUE 

945  F0RHAT(/»I10f3E14.6) 

C Post  FRONT  - Reset  ell  eleeenet  connection  nueber  to  positive 

C vilues  for  subseouent  use  in  stress  calculation* 

DO  520  I^liNELEH 
DO  520  J*l>9 

520  LNODSdf  J)>IABS(LNODS(I>J)) 

DO  30  I>>2»13»2 
J«=I/2 

QVl(I-l)=ASDIS(LN0DS(ISUPRlfJ)»2-l) 

QVKI)  «ASDIS(LN0DS(ISUPR1>J)«2) 

QV ( I - 1 ) =ASDI S ( LNODS ( I SUPR2 i J ) *2- 1) 

30  QV(I)  =ASDIS(LN0DS(ISUPR2»J)»2) 

DO  32  K»lf4 
32  SIF(K)=0.0D0 

C Calculate  the  stress  intensity  factor  for  Super-Eleaent  nueber  *1* 

DO  36  I-l>18 

SIF(1)=SIF(1)+0V1(I)»BCRT(1»1»I) 

SIF(2)=SIF(2)+QVl{I)»BCRT(l»2fI) 

SIF(3)=SIF(3)+0V(I)  »BCRT(2flfI) 

36  SIF(4)*SIF(4)+QV{I)  »BCRT(2i2.I) 

WRITE(6.922)ISUPR1 
WRITE(6f923)  SIF(l)  » SIF(2) 

WRITE(6>922)ISUPR2 
WRITE(6»923)  SIF(3)  » SIF(4) 

922  F0RMAT(///f23X» 'STRESS  INTENSITY  FACTOR  FOR  SUPER- ELEMENT  #'iI3»/) 

923  F0RMAT(/.12X»'S.I*F.  (OPENING  NODE)-' #E12.4» 10X» 

♦ 'S.I.F.  (SHEAR  H0DE)='»E12.4f///) 

RETURN 

END 


APPENDIX  F 


STRESS- CONCENTRATION  COMPUTER  PROGRAM 

The  following  computer  program  is  based  on  Savin's 
closed- form  analytical  solution  of  Chapter  2.  It  is  used 
to  obtain  all  tables  and  figures  in  the  first  section  of 
Chapter  5. 

The  subroutine  (POLRT)  which  is  used  to  calculate 
the  roots  of  the  characteristic  equation  of  Chapter  2, 
Equation  (2.11),  is  not  listed  here  since  it  is  already 
listed  in  one  of  the  finite-element  computer  subroutines 
in  Appendix  E.  The  only  modification  of  the  subroutine 
POLRT  needed  is  to  convert  it  from  single  precision  to 
double  precision  in  order  to  be  applicable  to  the  follow- 
ing computer  program. 
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IMPLICIT  REAL*8(A-H»0-Z) 

DIMENSION  XC0F<5) » SF ( 6 ) » COF ( 5 ) » ROOTR ( 4 ) » ROOTI < 4 ) 

REAL*4  SS(183) »XX(183) 

CALL  PLOTS( 10.0» 12.0»0f 0.5» 1 .5) 

XX( 182)=0.0 
XX( 133)=20.0 
SS< 132)=-3.0 
SS(183>=  10.0 
CALL  LINEWT(-l) 

CALL  AXIS(0.0»0.0» 'ANGLE' » -5 » 9 . 0 > 0 . 0 » XX ( 1 82 ) » XX ( 1 83 ) » 20 . 0 ) 
CALL  AXIS(0.0»0.0> 'STRESS-CONCENTRATION  FACT0RS'>28> 
*8.0»90.0fSS(182) »SS( 183) » 20.0) 

CALL  PL0T<9.0f 0.0»3) 

CALL  PL0T(9.0»6.0.2) 

CALL  PL0T(0.0»6.0»2) 

CALL  SYM<4.6»-0.4» . 1 7 » 280 > 0 . 0 » - 1 ) 

CALL  SYM( -0 .3>4.6»0.17»'l4SI7l4YI8l4/llP'f90.0»16) 

CALL  SYM<0.3»5.5f 0. 15f ' I OCJ I 1 ISOTROPIC' »0.0» 15) 

CALL  SYM(0.3»5.2»0.17» ' 10  J I 30  = 0 I 0= ' f 0 . 0 » 1 1 ) 

CALL  SYM(0.3»4.9»0. 17» ' I OA J I 30  = 30  I 0= ' » 0 . 0 » 1 3 ) 

CALL  SYM(0.3»4.6»0. 17f ' I OB J I 30  = 45  I 0= ' »0.0»13) 

CALL  SYM(0.3»4.3»0.17»'I0EJI30=60I0='»0.0»13) 

CALL  SYM(0.3»4.0»0. 17» ' I ODJ I 30  = 90  I 0= ' »0.0»13) 

DO  95  K=l»6 

M = 4 

IF<K.E0.6)  GO  TO  18 
RATI0=9.0D0 

READ(5»1005)  YOUNGl » Y0UNG2 f PO I SI  2 f G 1 2 » ANGL 

WRITE (6 f 1006) YOUNG 1 » Y0UNG2 » PO I S 1 2 » G1 2 » ANGL 

S11=1/Y0UNG1 

S12=-P0IS12*S11 

S22=1/Y0UNG2 

S66=l/G12 

* 

ANGL=ANGL*DATAN( 1 .0D0)/45 
SIN1=DSIN(ANGL) 
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C0S1=DC0S(ANGL) 

SIN2=SIN1*SIN1 

COS2=COS1*COS1 

SIN4=SIN2*SIN2 

C0S4=C0S2*C0S2 

3IN3=SIN2*SIN1 

C0S3=C0S2*C0S1 

3P( 1 )=S1 1*C0S4+(2*S12+S66)*SIN2»C0S2+S22*SIN4 
3F(2)=312*(3IN4  + C034)  + (311+322-366)!(t3IN2*C032 
3F'(3)=S11*3IN4+  (2*312  + 366)*3IN2*C032  + 322*C0S4 
3P(4)  = ( 2*311 -2*3 12-366)*3INl*C033-(  2*322-2*3 12-S66)*S I N33CC031 
3P(5)  = ( 2*3 11-2*31 2-366) *3 1 N3*C031-( 2*322-2*31 2-366) *3IN1*C033 
3P(6)=2* (2*311 +2*322-4*312-366 )*3IN2*C0S2+366*(3IN4+C034) 
URITE(6»9) (3P( I ) f 1=1 »6) 

8 FORMAT (2X» ' SP ( I ) = ' » 6E 1 4 . 6 f // ) 

XCQF( 1 )=SP(3) 

XC0F(2 )=-2*3P(5) 

XCQF(3)=2*3P(2)+3P(6) 

XC0F(4)=-2*3P(4) 

XC0F(5)=SP( 1 ) 

1C05  F0RMAT(5D12 .0) 

1006  FORMAT (/»3Xf ' YOUNG  1 = ' , F 1 2 . 3 » 8X » ' Y0UNG2= ' » F 1 2 . 3 » 8X » ' PO 1 3 1 2= ' r F7 . 5, 
» 6X»  'G12=' »F12.3»6Xf 'ANGLE=' »F6. 1 »//) 

CALL  POLRT(XCOF.COF»M»ROOTRf ROOT! f lER) 

DO  7 1=1,4 

7 WRITE(6,10)  ROOTR( I ) ,ROOTI ( I ) 

10  F0RMAT(5XfF10.5f 6X,F10.5, 'i ' ,/) 

U=R00TR(2)*R00TR(4)-R00TI (2)*R00TI (4 ) 
U=ROOTR(2)*ROOTI(4)+ROOTR(4)*ROOTI<2) 

X=R00TR(2)+R00TR(4) 

Y=R00TI (2)+R00TI (4) 

IF(K.LT.6)  GO  TO  19 
U=-l .ODO 

v=o.ono 

X=O.ODO 
Y=2.0D0 
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19  U2=U*U 
V2=V*V 
X2=X*X 
Y2=Y*Y 

Zl=2-U-X2-Y2 
Z2=2»X*U-2*X+2*Y*V 
Z3=1+2*U-U2-V2 
DO  20  I=lflSl 

THETA=(I-1)*I'ATAN(1.0D0)/45 

IF( ( I-l ) .EQ. 130)  GO  TO  25 
A=DSIN(THETA) 

GO  TO  26 

25  A=O.ODO 

26  A2=A*A 
A3=A*A2 
A4=A*A3 
A5=A*A4 

IF< ( I-l ) .EQ.90)  GO  TO  30 
BP=DCOS (THETA ) 

GO  TO  31 

30  EF'  = O.ODO 

31  B=RATIO*EF 
B2=E*B 
B3=B^B2 
B4=B*B3 
B5=B*B4 
Z4=A2*Y-A*B*V 
Z5=A2-A*B*X+B2*U 
Z52=Z5*Z5 
Z6=(B*V-A»Y)*B 
Z62=Z6*Z6 
Z7=A2+B2 
Z72=Z7*Z7 
Z8=A*B4*X-B5*U 

ZZ=A5*X+A4*B*Z1+A3*B2*Z2+A2*B3*Z3-Z3 
STR=< A2+RATI0*(BP*ZZ+Z4*Z72)/(Z52+Z62) )/Z7 
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XX(I)=I-1 
SS< I )=STR 
20  CONTINUE 

DO  97  I=lrl85»10 

11=1+9 

IFd.GT.ieO)  11  = 183 

97  WRITE (6»96) (XX( IN) » IN=I 
DO  98  I=lfl85»10 
11=1+9 

IF(I.GT.180)  11=183 

98  WRITE(6»9A) (SS< IN) f IN=I 
96  F0RMAT(3Xf 10F10.3f /) 

KKK=K-1 

CALL  LINE(XX»3S» 181 »1  »5 
95  CONTINUE 

CALL  PL0T(0.0»0.0»999) 

STOP 

END 


» II  ) 


» 1 1 ) 
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